Solving the nonlinear Schrödinger equation using exponential integrators

This PhD-thesis contains an introduction and six research papers sorted chronologically, of which the first four are accepted for publication. The introduction aims at giving a very brief summary of the background theory needed for the following papers. Also, some motivation of the issues addressed by the papers is given. Paper I discusses algebraic structures of ordered rooted trees and their applications to Lie group integrators. Results from Hopf algebra theory on elementary differentials for Lie group integrators are used, and applications to order analysis and backward error analysis are given. Paper II, III, IV, and V are primarily on exponential integrators, a class of numerical schemes tailored the solution of stiff systems of systems of ordinary differential equations. Paper II discusses classical order analysis and gives some theoretical results on the form of the integrators, applicable for the construction of high order exponential integrators. Paper III is on an implementation of exponential integrators on computers, and source code, available electronically, accompanies the paper. Paper IV includes an analytical and numerical study of the performance of two classes of exponential integrators on the nonlinear Schrodinger equation. Paper V is a numerical study of behaviour over long integration invervals on the nonlinear Schrodinger equation, using nonlinear spectral theory for determining validity of the numerical solution and thereby jugdging the numerical integrators. At last, in Paper VI, properties of a class of exponential like functions, essential in exponential integrators, are derived, using an approach based on Lie group theory.


Introduction
Although not new, exponential integrators were not considered a practical means of resolving systems of ordinary differential equations until very recently.Exponential integrators are especially designed to handle stiff systems, and accomplish this goal by constructing exact integral curves for the linear part of the differential operator.Constructing the integral curves entails the application of the matrix exponential and related functions.
The class of integrators henceforth termed exponential integrators first appeared in Certaine [5] and Nørsett [16].These schemes are both of exponential time differencing (ETD) type.Then Lawson [14] constructed the integrating factor type in 1967.Recent reports on exponential integrators show that especially for parabolic semi-linear problems, the ETD type of exponential integrators outperform integrators of Lawson type [13,15,17].However, few results are avail-able with respect to the performance on non-parabolic problems like the Schrödinger equation.
In this paper we test a fourth order Lawson integrating factor scheme against a fourth order ETD scheme, ETD4RK in [6].Most other similar exponential integrator schemes perform very similarly to the ETD4RKscheme.Exponential integrators are introduced in Section 2 and some analysis and numerical results are presented in Sections 4 and 5.
The equation we will use for numerical tests in this paper is the nonlinear Schrödinger equation in one space dimension for all t ≥ 0 ψ(x, 0) = ψ 0 (x), x ∈ [−π, π]. ( This Schrödinger equation arises in several different areas of physics of which we mention multiscale perturbation theory, gravity waves in water, and propagation of intense optical pulses in fibres.The nonlinearity constant λ controls the ratio of dispersive effects to nonlinear effects, and may give a focusing version of the equation.The equation may be both parabolic and hyperbolic, it has some smoothing effects, but timereversibility prevents it from generating an analytic semigroup, which is fundamental for the stiff order analysis in Section 3.An introduction to the mathematical theory of the nonlinear Schrödinger equation is given in [4].
We would like to point out that we do not try to directly preserve any invariants of the equations in question, as opposed to many other specialized schemes for the Schrödinger equation.In this work, we test the given schemes on a limited time scale, and focus on reporting the observed order.The Schrödinger equation possesses several conservation laws, notably conservation of density, energy and momentum.For long-time integration where stability and preservation of invariants is an important factor, multisymplectic schemes may be a viable choice [3,10].The benefits of preser-vation of invariants must be weighted against the additional cost necessary for multisymplectic schemes.
For our Schrödinger equation (1) we will employ a discrete Fourier transform with N F modes.Upon semi-discretizing the physical problem in space, we obtain a system of ordinary differential equations given by ẏ = Ly + N(y) in which y ∈ C n is the Fourier transform of ψ, L ∈ C n×n , and N : C n → C n .For the Schrödinger equation (1) the L matrix becomes diagonal with entries and the nonlinear function N(y) becomes in which each component of y represents a particular Fourier mode, k.

Exponential integrators
Exponential integrators are explicit schemes which recover the exact solution to linear problems.As such, this class of schemes is well suited to problems which can be split into a linear and a nonlinear part, and for which the linear part is either stiff or unbounded and the nonlinear part grows more slowly than the linear part.When semi-discretizing PDEs, this happens if spatial derivatives in the linear part are of higher order than in the nonlinear part.We note that the Schrödinger equation, whether semi-discretized as in (2) or in its original PDE form (1), satisfies these requirements although the matrix L of ( 2) is unbounded only in terms of the parameter N F .
In the following, we consider systems of ordinary differential equations split into a linear and a nonlinear part as ẏ = Ly + N(y,t), y(0) = y 0 . ( As alluded to in the above paragraph, exponential integrators applied to this problem possess two primary features 1.If L = 0, the integrator reduces to a classical Runge-Kutta or linear multistep method. 2. If N(y,t) = 0 for all y and t, the integrator reproduces the exact solution to (5).
The nonlinear function N may depend on time, but the linear part should not be explicitly time dependent in order for the exponential integrator to be computationally competitive.Moreover, exponential integrators implicitly assume that most of the system's inherent dynamic behaviour can be ascribed to the linear operator L.
Classical integrators are divided into two classes; linear multistep methods and one-step Runge-Kutta methods.This paper considers only exponential Runge-Kutta methods.We note that the framework of general linear methods, a generalization of both linear multistep methods and Runge-Kutta methods, may also be extended to define exponential integrators as in [2].
Exponential integrators of Runge-Kutta type are written as in which Y i , i = 1, . . ., s are internal stages and y 1 is the final approximation of y(t 1 ) = y(t 0 + h).This format extends the common format of Runge-Kutta schemes in that the coefficients a i j and b i are now analytic functions of the linear operator L.
In order to fulfill the two features of an exponential integrator, a i j (0) and b i (0) must be the coefficients of some underlying Runge-Kutta-method.It is evident that this scheme will solve linear equations (N(y,t) = 0) exactly.Extending the notation of Butcher, the coefficient functions and collocation nodes are written up in the tableau where we have used z = hL for convencience.The two simplest choices of exponential integrators of Runge-Kutta type are the Lawson-Euler 0 0 e z and Nørsett-Euler schemes.The function ϕ 1 (z) in the Nørsett-Euler scheme is given by ϕ 1 (z) = (e z − 1)/z.The latter scheme has been reinvented several times, and is also known as ETD Euler, filtered Euler, Lie-Euler (using the affine Lie group) and exponentially fitted Euler.

Lawson schemes
The Lawson exponential integrators, of which Lawson-Euler is a special case, are derived by introducing a change of variables involving an integrating factor and applying a classical Runge-Kutta scheme to the transformed equation.Given an underlying Runge-Kutta scheme with coefficients ãi j , bi and corresponding quadrature nodes c i , the Lawson exponential integrator coefficient functions are as given in [14] Lawson schemes are particularly simple to implement, but have some disadvantages as reported early in the history of exponential integrators.For example, they do not preserve fixed points of the differential equation, and are also known for rather large error constants.
The aim of this paper is to elaborate on the performance of a Lawson exponential integrator based on Kutta's classical fourth order method.This scheme is given by the tableau and will be denoted "Lawson4".Ehle and Lawson modified the Lawson schemes in their paper [7] and introduced another fourth order exponential integrator also using the ϕ 1function, thereby slightly improving the performance for parabolic applications and regaining fixed point preservation.Their modification was in the direction of ETD-schemes, but it is not competitive to the now known ETD-schemes.

Exponential time differencing (ETD)
Rather than using integrating factors, we may approximate the nonlinear function N(y,t) by some polynomial in t, and integrate the approximate equation exactly.The resulting schemes are known in recent literature as the "exponential time differencing" (ETD) schemes, although this name is not entirely descriptive.The polynomial approximation may be calculated using previous steps of the integration process, thus producing multistep ETD schemes, or by Runge-Kutta-like stages, resulting in ETD schemes of Runge-Kutta type.We refer the reader to the review paper [15] for a thorough review of exponential integrators of these types.
For notational simplicity, and without loss of generality, we consider only autonomous problems N(y) = N(y(t)) in the remainder of this paper.
Lemma 2.1.The exact solution of the initial value problem has the expansion where Proof.The basic idea is just a Taylor expansion of the nonlinear function N(y(t)) and the variation of constants formula.A proof may be found in [13, Lemma 1.1].
We will in this paper compare the Lawson4 scheme (11) against the most commonly used fourth order ETD scheme, ETD4RK, due to Cox and Matthews [6].The coefficients of ETD4RK are given by 0 in which Computationally, the Lawson4 scheme ( 11) is much cheaper and easier to implement on a computer than ETD4RK.The evaluation of ϕ-functions in (12) has numerical issues, and we believe this is best dealt with using scaling and corrected squaring together with Padé approximants.Details on this may be found in [2].
3 Order conditions Classical order analysis for numerical integrators develops Taylor expansions for all quantities.The analysis, however, is rigorous and valid only in the limit as hL → 0. If L is defined by spatially semi-discretizing an unbounded differential operator L , L may be unbounded in terms of a parameter, typically the spatial resolution.Thus, hL → 0 cannot generally be guaranteed independently of the parameter.As such, classical order analysis is of somewhat limited use in the study of exponential integrators applied to unbounded semilinear problems.Nevertheless, classical order conditions must be satisfied for exponential integrators and traditional Runge-Kutta integrators alike.The Law-son4 (11) and ETD4RK (13) schemes are methods with classical order four.Details on classical order analysis for exponential integrators using B-series may be found in the paper [1].
A recent paper of Hochbruck and Ostermann [9] studies exponential integrators applied to infinite dimensional semi-linear parabolic Cauchy problems.Conditions under which the integrators converge in this abstract setting are rather restrictive, and give rise to the notion of stiff order.Including the classical order conditions as special cases, these "stiff order conditions" constitute an extended set of requirements which must be satisfied to guarantee high convergence rates.In this context Lawson4 (11) has stiff order only 1 and ETD4RK (13) has stiff order only 2. The use of ϕ-functions in the coefficient functions a i j (z) and b i (z) of ( 6) is required to attain high stiff order.
However, the applicability of stiff order analysis to the nonlinear Schrödinger equation remains an open issue.Integrators of stiff order four, examples of which are listed in [2], perform similarly to ETD4RK in this study.This suggests that high stiff order is not critical to achieving efficient schemes in all cases, and these high stiff order schemes are therefore omitted in all plots.
The first stiff order condition for an exponential integrator is easily obtained by comparing the numerical solution given in (6) to the exact solution from Lemma 2.1.For the first order in h we get the equation which we rewrite as Based on this, the first stiff order conditions reads The Lawson integrators do not satisfy this condition exactly, but the integrators nevertheless satisfy the condition to a sufficient degree of accuracy, a notion which will be explained in Section 4.There we study the solution's dependence upon the Schrödinger equation potential function V (x).
An easy route to deriving two stiff order conditions is considering preservation of fixed points.Exact preservation of fixed points is important in many applications, and hence a desirable property of exponential integrators.Requiring Ly = −N(y) and y 1 = y 0 , equation (6b) gives For equation (6a) we require Y i = y 0 for all i, and we obtain These are precisely the first and third stiff order conditions in [9].Lawson integrators fulfill neither of these conditions, and thus do not preserve fixed points.ETD4RK, however, fulfills both conditions for fixed point preservation.
Despite their low stiff order, Lawson4 (stiff order 1) and ETD4RK (stiff order 2) still behave as fourth order schemes on our problem, given smooth initial condition and smooth potential.See Figure 1.

Potential function dependency
The first stiff order condition (15) is not satisfied by the Lawson schemes.The significance of the stiff order conditions in the case of non-parabolic problems like the Schrödinger equation is unclear, but the conditions still affect the numerical performance in some cases.Figure 1 shows that the Lawson scheme is rougly 100 times more accurate than ETD4RK at comparable stepsizes, however, as we will justify, the performance results in Figure 1 are strongly influenced by the smoothness of the potential function used in this particular test.
In this section we study how the regularity of the potential function V (x) affects the numerical performance of the Lawson4 integrator.

Order estimates
The analysis will be based on the rate of decay for the Fourier coefficients of input functions.The relationship between Fourier decay and differentiability is taken from a well-known result in Fourier analysis.
We estimate the error contribution from the first stiff order condition (15) in Fourier space.By substitution of the b i (z) of Lawson4 into (15) we obtain which, when z is small, has the Taylor expansion For x ∈ R we have |ψ 1 (xi)| ≤ 2, and we use this to construct an upper bound for ψ 1 for high Fourier modes where the Taylor expansion ( 19) is not valid.Let the bound be the function which is sufficiently sharp for our purpose.2) has a decay rate of at least r for all time and 1 2 < r ≤ 2p, that is

N(y(t)) (in Fourier space) in (
then the local error contribution for the first stiff order condition is Proof.We bound the ψ 1 (−ik 2 ) function by where the critical mode value is k c = (2C) To estimate the error, we sum over k in the first stiff order condition in which we estimate the sums using Euler-Mac-Laurin's summation formula we get the dependency on h, and the square root of the dominating term is If r > 2p, the scheme is not accurate enough to capture the "non-smoothness" of the N-function, and the first order condition does not contribute to any error of order less than the classical error.
Thus, as long as the nonlinear function is smooth enough, we can also include Lawson4 as one of the schemes that obey the first stiff order condition, although only accurately enough so that its main features as a fourth order classical method is conserved.Looking at only the first stiff order condition is sufficient for explaining the observed numerical behaviour in this paper.

Numerical results
In the following experiments, we have used an artificially constructed potential with a prescribed decay rate r.This means constructing the potential by letting its Fourier modes be 1/(ik r ) multiplied with a complex number in which both the real and the imaginary part are normally distributed with mean zero and variance one.Then we have used MATLAB's inverse discrete fourier transform to get an example function for use.We note that in particular the hat function has decay rate of 2, although Lemma 4.1 only predicts 1.This is due to bounded variation of the hat function.
Figures 3 and 4 show observed error behaviour when solving the nonlinear Schrödinger equation subject to a smooth initial condition and potential functions of regularity 2 and 4 respectively.Low regularity   potential functions lower the regularity of the nonlinear function N(y(x,t)).Assuming N is no more regular than V (x), Proposition 4.2 then predicts orders 1.75 and 2.75 respectively for the Lawson4 scheme in these cases.We conclude that the observed order corresponds fairly well to what is predicted by the proposition.Moreover, we see from the plots that for the Lawson schemes, the global error as a function of time step oscillates rather wildly when not all eigenmodes are resolved by a small enough h.These oscillations are smooth on a zoomed plot and are due to some resonance effect.This is further discussed in Section 6.2.
In this section we will see that the ETD4RK scheme is more influenced by the regularity of the initial condition than is the Lawson4 scheme.A crucial introductory numerical observation is that the dependency on the initial condition is present in the linearized version as well as in the nonlinear version, that is when λ = 0 in (1).This facilitates substantially simplified analysis.where L is the Laplacian in Fourier domain as before (diagonal, ik 2 ) and W is a circulant matrix stemming from a Fourier transform of the potential in the Schrödinger equation (1).The fact that the matrices L and W in general do not commute is the source of the order reduction observed in the Lawson scheme as we shall see.If, on the other hand, the potential function is a constant, L and W will commute, and order reduction is not observed.

Analysis for the linear problem
The presentation here resembles the Strang splitting analysis of Jahnke and Lubich [11] on (22) when the linear operator in the differential equation is unbounded.They found order reduction due to the same phenomena that we will see here.
Applying an explicit exponential integrator to (22) we get Then, by way of the variation of constants formula, the exact solution to (22) may be represented as e h(L+W ) y 0 = e hL y 0 + h 0 e sL W e (h−s)(L+W ) y 0 ds.
For our fourth order schemes, we iterate the variation of constants formula four times for the exact solution resulting in a sum including up to fivedimensional integrals.Applying the variation of constants formula once more to remove W from the exponential, and substituting θ = (h − s)/h, it is clear that a second order scheme must satisfy

Regularity requirement for the Lawson scheme
Inserting the Lawson4 scheme b i (hL) coefficients allows an immediate interpretation of (24).The left hand side of ( 24) is Simpson's quadrature of the function f (θ ) = e (1−θ )hL W e θ hL .The error of Simpson's quadrature is known to be f (4) (ξ )/2880 for some ξ ∈ [0, 1], and in this case Transforming from Fourier space to phase space, L becomes d 2 /dx 2 and W becomes a multiplication operator denoted by V .One may verify the formula When m = 4, one observes that the Lawson4 scheme satisfies condition (24) to a sufficient degree of accuracy if the initial condition in phase space ψ 0 (x) ∈ C 4 (−π, π) and the potential V (x) ∈ C 8 (−π, π).
Iterating the variation of constants formula further, one obtains additional iterated integrals.As these integrals involve only lower derivatives of the appropriate f (θ 1 , θ 2 , . . . ) function, equating to lowered regularity requirements for V and y 0 , we omit the details in this exposition.

Regularity requirement for ETD4RK
We interpret (24) in a Gauss quadrature sense with the weight function w(θ ) = e (1−θ )hL .Requring the quadrature formula to be exact for fourth degree polynomials gives four stiff order conditions.
For ETD4RK this is not in general satisfied when k = 4, and we expect the principal quadrature error term to depend on g (4) (θ ) where g(θ ) = W e θ hL u 0 .Differentiating this function, we get g (4) (θ ) = h 4 W L 4 e θ hL y 0 (28) an upper bound for which translates to y 0 being at least 8 times continuously differnentiable in space.Thus, we should expect ETD4RK to demand more regularity for the initial condition than Lawson4.On the other hand, ETD4RK makes no demand on the regularity of the potential function, as opposed to the Lawson4 scheme.
Figure 7: Lawson4 is close to order 4 but oscillating in the linear case, λ = 0, while ETD4RK suffers from the low regularity (2) of the initial function.The potential is smooth, 1/(1 + sin 2 (x)).In the corresponding local error plot, ETD4RK exhibits order 1.75 and Lawson order 5 with no oscillations.
ETD type schemes.Particularly, the coefficient functions of Lawson schemes are given explicitly by (10), whereas derivation of ETD type coefficient functions is typically more cumbersome.Additionally, Lawson type schemes require one or more matrix exponentials for which acceptable algorithms are well known.
The ETD type schemes require the evaluation of multiple ϕ-functions, a computational problem which is at least as difficult as computing matrix exponentials.
In evaluating ϕ-functions, Kassam and Trefethen [12] discovered a stability problem which they solved by contour integral evaluation in the complex plane.This requires an a priori contour radius which in general is problem dependent and not trivially available.In our numerical experiments, we found a scaling and squaring technique together with Padé-approximations of the ϕ-function to be a better option, inspired by a code from [8].The actual implementation is discussed in [2].

Oscillations in observed order
Most order plots for the Lawson4 integrator show significant oscillations in observed accuracy as a function of timestep h.Zooming in on each plot reveals that the the oscillations are smooth but quickly varying magnitudes of the highest eigenmode of L. These oscillations span roughly 2 orders of magnitude, and therefore represent a considerable error contribution at particular time step sizes.
The oscillations are due to some resonance effects, and that these are not damped as in the case of ETD schemes by dividing by z in the ϕ-functions.To avoid these oscillations the Lawson schemes therefore must use ϕ-based coefficient functions.This, in turn, effectively renders the scheme into another type than what has been denoted Lawson schemes in this paper.Moreover, in a sense the resulting scheme is worse than Lawson's scheme as the modified scheme becomes more sensitive to the regularity of the initial condition.

Low regularity potential and initial condition
Using low regularity initial conditions and potential functions, we get the mixed case of undesirable behaviour from both types of schemes.Varying both the initial condition and the potential (one particular combination of which is shown in Figure 8), there is little actual gain from choosing one scheme over the other.However, due to the observed oscillations, ETD4RK might be a better choice in these cases.

Exponential general linear methods
General linear methods generalize Runge-Kutta integrators and multistep integrators.Exponential general linear methods thus generalize the integrators catered for in this paper, as reported in [2].A class of exponential general linear methods known as the generalized Lawson schemes, see [17], mixes the Lawson and ETD schemes and give good results on parabolic problems, achieving high stiff order.However, in the experiments described in this paper, these schemes never perform better than the best of ETD4RK or Lawson4.

Rounding error accumulation
In closing we would like to comment on an important feature of our experiments.The measured error does not decrease further as a function of decreasing stepsize once the error reaches a level of 10 −10 .As this is several orders of magnitude larger than machine accuracy, it is clear that rounding errors introduced in the evaluation of the ϕ-functions affect long-time accuracy of the exponential integrator.We still believe that the Padé approximation, as described in [2], is the best algorithm for evaluating ϕ-functions, and that accuracy of exponential integrators may be increased by further research into this algorithm.

Conclusion
We have studied the numerical performance of the Lawson4 scheme compared to the ETD4RK scheme on a nonlinear Schrödinger test problem and observe that the actual performance is heavily influenced by the potential function and initial condition.In short, Law-son4 is dependent upon the regularity of the potential function while ETD4RK is dependent upon the regularity of the initial condition.Stiff order conditions are used as a tool for explaining the observed behaviour, although the general applicability of stiff order conditions to non-parabolic problems remains unclear.Further research is necessary to explain phenomena exhibited by exponential integrators on partial differential equations.

Figure 1 :
Figure 1: A global order test.Both exponential tegrators in this study perform as order 4 integrators.The dotted line is only an indicator line showing how order 4 looks like.

Figure 2 :
Figure 2: The error in the first stiff order condition for Lawson4, h = 0.1 in this plot.

Figure 3 :
Figure 3: Global error when the potential has regularity 2.

Figure 4 :
Figure 4: Global error when the potential has regularity 4.

Figure 8 :
Figure 8: The mixed case: Both the initial condition and the potential has regularity 2.