Research Article  Open Access
Large TimeStepping Spectral Methods for the Semiclassical Limit of the Defocusing Nonlinear Schrödinger Equation
Abstract
We analyze a class of large timestepping Fourier spectral methods for the semiclassical limit of the defocusing Nonlinear Schrödinger equation and provide highly stable methods which allow much larger time step than for a standard implicitexplicit approach. An extra term, which is consistent with the order of the time discretization, is added to stabilize the numerical schemes. Meanwhile, the firstorder and secondorder semiimplicit schemes are constructed and analyzed. Finally the numerical experiments are performed to demonstrate the effectiveness of the large timestepping approaches.
1. Introduction
Consider the initialvalue problem for the semiclassical limit of defocusing nonlinear Schrödinger equation (DFNLS) in one space dimension where , the domain is an open set in , is the complexvalued wave function, is the scaled Planck constant, and is a real constant. is given with the initial amplitude and is the real initial phase function with and decaying rapidly as . To avoid boundary effects and obtain enough information of statistics, we compute the solution in a sufficiently large interval .
The DFNLS (1.1) has infinitely many conservation laws (see, e.g., [1–3]), two of which are
where (the position density) and (the current density) are primary physical quantities and can be computed from the wave function :
where ““ denotes complex conjugation. The two conservation laws can also be denoted by the norm and energy density conservation laws
where is the standard norm in .
The DFNLS is connected with many applications in science and technology. One of the most important applications of the semiclassical limit of (1.1) arises in nonlinear optics [4], in particular regarding the production and propagation of laser pulses in fiber optics. Here depends on the dimensions on the problem (amplitude of the initial pulse, dimensions of the fiber, etc.). Another application appears in the study of the socalled BoseEinstein condensate [5, 6], where now does represent the Planck constant.
The limit is called the semiclassical limit and considerable attention has been given recently to the investigation of its existence and structure [7–11]. The dynamics of the limit is an open problem. Wright et al. have provided the exact solution of the geometric optics approximation of the defocusing nonlinear Schrödinger equation in [12]. Much progress has been made recently in understanding semiclassical limits of the linear Schrödinger equation. Particularly by the introduction of tools from microlocal analysis, such as defect measures [13], Hmeasures [14], and Wigner measures [15, 16].
Numerically this is a notoriously difficult problem. The oscillatory nature of the solutions of the nonlinear Schrödinger equation with small provides severe numerical burdens. Even for stable numerical approximations (or under mesh size and time step restrictions which guarantee stability), the oscillations may very well pollute the solution in such a way that the quadratic macroscopic quantities and other physical observables come out completely wrong unless the spatialtemporal oscillations are fully resolved numerically. For a longtime considerable consideration has been given to this problem by many researchers. In [17, 18], Markowich et al. utilized the Wigner measure in analyzing the semiclassical limit for the the linear Schrödinger equation with small , to study the finite difference approximation. In [19, 20], Bao et al. simulated the solution by the timesplitting spectral approximation for the linear and nonlinear Schrödinger equation. Hector D. Ceniceros presented the semiclassical limit of the focusing nonlinear Schrödinger equation by moving mesh method [21] and Ceniceros and Tian studied it by nonlinear Fourier filtering method in [22]. But there are few research on DFNLS.
The purpose of this paper is to provide an effective numerical method to solve the DFNLS with small parameter . In space direction, we use Fourier spectral method, and in time direction we use difference method. In the case that the space direction has spectral accuracy, we expect to improve the computational efficiency by increasing the time step, which is beneficial to the numerical simulation and the discussion on the longtime behavior of this problem. In numerical computation, for one thing, because of the strong nonlinear character of (1.1), the numerical stability cannot be completely ensured by the traditional numerical methods. For another thing, because of the oscillation of high frequency of the solution of this equation, an ideal discretization of space and time should decompose bordering surface precisely. Hence a high resolution ratio of space is needed, and very small time step must be used to satisfy the demands of computational stability if we take standard explicit or semiimplicit scheme. When is very small, the explicit treatment for nonlinear item usually causes strict restriction on time step, and it is not practical to have a longtime computation in this case. To relax restriction of this kind of problem, we provide a semiimplicit difference Fourier spectral method. That is to say, to improve the stability of algorithm, we add a routine item to the semiimplicit scheme where a larger time step is allowed in the numerical experiment [23–25]. First, for the treatment of time discretization, the semiimplicit scheme is used. Second, for the secondorder item in (1.1), the implicit scheme is employed to reduce the restriction on time step. Third, for the nonlinear item, the explicit scheme is used to avoid solving systems of nonlinear equations at each time step as this explicit treatment can be computed by FFT. We expect that the computation will be easy, large time step is allowed in the solution procedure, and the numerical stability can be guaranteed at the same time.
The organization of this paper is as follows. In Section 2, a theoretical analysis for the semiimplicit method is provided, and the firstorder timestepping method stability properties are investigated. The secondorder semiimplicit methods are mainly investigated in Section 3. Full discretization schemes and accuracy tests are presented in Section 4. Finally, some computational results are given in Section 5.
2. Stability Analysis for FirstOrder Time Discretization for (1.1)
In this section, we present Fourier spectral approximations of the problem (1.1) with periodic boundary conditions.
We define periodic Sobolev space as
Define inner product and norm as
For any let where
First we consider the classical firstorder semiimplicit method, the scheme is obtained by combining firstorder backward difference (BD1) for and a firstorder extrapolation (EP1) for the nonlinear item at the same time where is the timestep and , is an approximation to .
We expect that the implicit treatment for the secondorder item in (2.4) can relax the restriction on time step. And we notice that it will not increase the complexity of the algorithm through the spectral method because the secondorder item here is complete linear. In Section 4, it is showed that large step is not allowed in scheme (2.4) when is small. That is to say, to guarantee the stability, a very small time step must be used to get the precise solution. To overcome this shortage, we add an item to the scheme and have the following modified scheme (BD1/EP1):
where is an undetermined positive constant, which will be given in the following results. The purpose of adding an item is to improve the stability of the original scheme, and a larger time step can be allowed in the computation.
As we can see, (2.5) can be denoted as the following weak form:
Theorem 2.1. If in (2.5) is sufficiently large, the following energy inequality holds: where More precisely, the positive constant satisfies
Proof. Let in (2.6), take the real part of the result
Now we estimate the terms in (2.10), respectively,
From the inequality
and the above estimations, we obtain
The last term mentioned previously can be made nonnegative if the following inequality holds:
Remark 2.2. From the proof of Theorem 2.1, we know that the selection of in (2.14) depends on solution of the equation in the th step. But essentially it is not difficult to choose a suitable . Actually, we can determine it approximately from the following expression:
Remark 2.3. In numerical computation, condition (2.14) can be guaranteed by a posterior method; that is, when computation of each step finishes, we check whether the condition (2.14) is satisfied to guarantee the decrement or the conservation of discrete energy.
3. HigherOrder Methods for SemiImplicit Time Discretization
3.1. SecondOrder Scheme: BD2/EP2
Similar to the firstorder scheme, HigherOrder time discretization can be constructed. For example, by combining the secondorder backward difference (BD2) for , and the secondorder extrapolation (EP2) for the explicit treatment of nonlinear item, we have the following secondorder scheme:
here we add an item similar to the item added in the firstorder scheme, and expect that this item can be consistent with the global order of the scheme and keep the scheme stable. Therefore an item is added in the above scheme, and we arrive at the following modified secondorder scheme (BD2/EP2):
As usual, to start the iteration, is given by the initial condition and is computed by the firstorder scheme (2.5).
Equation (3.2) can also be denoted by the following weak form: where is a positive constant.
Lemma 3.1. One has
Proof. We have then
Lemma 3.2. One has
Theorem 3.3. Let be the solution of the semidiscrete problem (3.2). The corresponding generalized energy is defined as follows: Then, for large enough , the following discrete energy inequality holds
Proof. Let in (3.3), taking the real part of the result we have
Using Lemmas 3.1 and 3.2, the first term of the left of (3.10) is
the third term of the left of (3.10) is
from
we have
Using the definition of and from the estimation of (3.11), (3.12), then
where
Obviously, (3.9) holds provided that . In the following we will show that can be made nonnegative if is chosen large enough. To simplify we denote by . By simple calculation and applying Young's inequality to the suitable term, we get
that is,
We choose such that the integral term in the above estimate is nonnegative, then constant holds
Under condition (3.20) on , we have
This completes the proof of this theorem.
3.2. ThirdOrder Scheme: BD3/EP3
thirdorder scheme for solving the semiclassical limit of the defocusing nonlinear Schrödinger equation can be constructed in a similar manner as used in the subsection. Specifically, we can obtain the BD3/EP3 scheme in the following form:
where, in order to start the iteration, are calculated via a first and secondorder scheme, respectively. The stability analysis of the scheme (3.22) requires some very detailed energy estimates and will not be presented here. The numerical results obtained in the next section indicate that the thirdorder time discretization of type (3.22) is also stable as long as the constant is sufficiently large. But comparing with the lowerorder scheme, the influence on the stability from item in the thirdorder scheme is not very notable.
4. Full Discretization Scheme and Accuracy Tests
A complete numerical algorithm also requires a discretization strategy in space. Since Fourier spectral method is one of the most suitable spatial approximation methods for periodic problems [26–28], here we get the space discretization of semiimplicit scheme from different time iteration schemes by this method.
Let
We use the following Fourier transformations:
The expansion of (2.5) is required to satisfy the following weak form:
Similarly, we have the full discretization Fourier spectral approximation of problem (3.2):
The space discetization of BD3/EP3 (3.22) has a similar form.
It deserves to mention that the corresponding spectral coefficient can be fully solved because all implicit items about in (4.3) and (4.4) are linear.
Choosing test function as the primary function of , we obtain a system of linear equations for each module in Fourier space:
The corresponding secondorder BD2/EP2 scheme is as follows:
Visibly solving (4.5) and (4.6) is efficient because we can solve the system of linear equations determined by (4.5) and (4.6) through obtaining the inverse of a simple diagonal matrix. For the problem of fully discretization (4.5) and (4.6) we derive two theorems similar to the energy inequalities in Theorems 2.1 and 3.3 (its proof will be omitted here).
Theorem 4.1. Consider the numerical scheme (4.3), if the solution of (4.3) satisfies where
Theorem 4.2. If then the solution of (4.4) satisfies where is positive constant and
Then we investigate the stability properties of different schemes by numerical experiments, where special attention will be paid to the improvements of stability caused by the item. In the computation we choose or as the domain of computation, and choose random number as the initial condition. First, we define as the largest time step which makes the computation stable; that is, when the time step is larger than the numerical solution will “blow up." For different , Table 1 lines up the values of in the Fourier spectral schemes of (2.5), (3.2), and (3.22), respectively, and the Fourier module used in the computation is . From observing the data in Table 1, we draw the following conclusion.

When is large, we see that the is also large. Since all standard semiimplicit schemes are stable, small item is not necessarily required. When is small, smaller time step is needed to satisfy the requirement of stability if we use a standard semiimplicit approach. This phenomenon is prominent especially for the case of higherorder or very small . For the firstorder scheme BD1/EP1, just as proved in Theorem 2.1, when is large enough the unconditional stability of the computational process can be guaranteed. For the secondorder scheme, item improves the stability obviously. Comparing with the case of , about 10 times timestep is allowed when we use item. For the thirdorder scheme, although the effect of adding item is not very obvious, adding the item scheme is still effective in improving the stability and enlarging the time step to a certain extent.
Now we turn to time accuracy comparison. Since the exact solution of (1.1) is unknown, we use numerical results of the BD2/EP2 with , and as the “exact" solution. We discuss it in the following two cases.
Case 1 (weak cubic defocusing nonlinearity). We set , the final time . Table 2 shows the errors obtained, at the same time we derive the rate of convergence about time by the following formulation: where is the error when using timestep , the rate may be verified experimentally by solving a problem on a sequence of finer and finer partitions and approximation. From Table 2, we can easily calculate , , , which is coincident with our theoretical conclusion.

Case 2 (strong cubic defocusing nonlinearity). Figure 1 shows the position density and the current density for , with different values of and in different time . It is observed that there is a good agreement between the numerical results obtained by using standard semiimplicit timestepping method (i.e., ) with small and the modified method (3.2) with larger .
(a)
(b)
(c)
(d)