WKB-method for the 1D Schrödinger equation in the semi-classical limit: enhanced phase treatment
This paper is concerned with the efficient numerical computation of solutions to the 1D stationary Schrödinger equation in the semiclassical limit in the highly oscillatory regime. A previous approach to this problem based on explicitly incorporating the leading terms of the WKB approximation is enhanced in two ways: first a refined error analysis for the method is presented for a not explicitly known WKB phase, and secondly the phase and its derivatives will be computed with spectral methods. The efficiency of the approach is illustrated for several examples.
Key words. Uniformly accurate scheme, Schrödinger equation, highly oscillating wave functions, higher order WKB-approximation, asymptotically correct finite difference scheme, spectral methods
AMS subject classifications. 35Q40, 81Q20, 65M70, 65L11
This paper is concerned with the numerical solution of highly oscillating differential equations of the type
where is a very small parameter and a sufficiently smooth function. Such models play an important role in electromagnetic and acoustic scattering (1D Maxwell and Helmholtz equations in the high frequency regime), wave propagation problems in quantum and plasma physics (Schrödinger equation in the semi-classical regime). For very small , the wave length is very small, such that the solution becomes highly oscillating. In a classical ODE–scheme (like in [9, 10]) such a situation requires a very fine mesh in order to accurately resolve the oscillations, see Fig. LABEL:onde. Hence, standard numerical methods would be very costly and inefficient here.
A possible remedy is to use analytic a-priori information on the solution to simplify its numerical treatment. In the semi-classical limit , asymptotically correct solution formulas are provided by the well-known WKB-approximation (cf. ). They can be obtained in arbitrary order of accuracy. For the numerical method to be discussed here, the second order WKB-approximation for (LABEL:EQ_ref) is relevant, and it reads
In recent years various numerical strategies for highly oscillatory problems from quantum mechanics have been developed. They are based either on adiabatic integrators (see [13, 12, 15]; §XIV of ), WKB-approximations (see [16, 1]), or a macroscopic reformulation . In this paper we shall extend and refine the asymptotically correct WKB-scheme from . Both of these papers are concerned only with the oscillatory case (i.e. ). The evanescent case (with ) is equally important in applications, but needs a different numerical approach, see . For the inclusion of a turning point we refer to .
In  the authors introduced a WKB-based method that allows for an efficient and accurate solution of (LABEL:EQ_ref) in the semi-classical regime. This is a hybrid method consisting of two steps: first an analytic preprocessing of (LABEL:EQ_ref) to transform it into a smoother, i.e. less oscillatory problem that can be solved numerically on a coarse grid with high accuracy. Then, the numerical solution of the transformed ODE problem is based on a truncated Picard iteration. This explicit representation of an approximate solution involves oscillatory integrals with two small parameters, and the step size . Hence, the key issue in the numerical step is to construct approximations of these integrals that are -uniform as well as accurate enough w.r.t. , since the latter determines the local discretization error for the ODE scheme. In  both a first and second order method (w.r.t. ) were derived.
In the analytic step, the second order differential equation (LABEL:EQ_ref) is first transformed to a system of first order differential equations using
The numerical analysis of the WKB-method in [1, Th. 3.1] led to the following error estimates for the first and, resp., second order methods:
where is a numerical approximation for the solution at the grid point . Here and in the sequel, denotes generic, but not necessarily equal constants that are independent of grid index , , and . Moreover, is the order of the chosen numerical integration method for computing the phase integral
which is a smooth function of both and . Here, denotes any vector norm in .
The second terms in (LABEL:error_U) and (LABEL:error_U_2ORD) are the errors due to the WKB method, and they are of the order , even for constant step size . By contrast, the first term is critical in the semi-classical limit as it grows for . This error is due to using an approximate phase for the analytic transformation step (back and forth). To limit the size of the first term one has two options: One can either choose an -dependent step size (like when using e.g. the Simpson rule with for the phase computation, cf. [15, 1]). Alternatively one can use a highly accurate method – like a spectral method, and this will be our approach here.
The goal of this paper is twofold: In , numerical errors in the phase were only considered for the backward transformation. Hence we shall first complete that error analysis by taking into account that also the first (analytic) transformation is typically affected by an inaccurate phase. Secondly, we shall combine the WKB-method with a spectral method for computing the phase function , yielding spectral accuracy. With little effort, this will reduce the first term in the error estimates (LABEL:error_U) and (LABEL:error_U_2ORD) to the order: machine precision over , making them irrelevant for most practical computations.
This paper is organized as follows: In §LABEL:SEC2 we briefly review the WKB-method from , and in §LABEL:SEC3 we extend the error analysis to include the phase error in . In §LABEL:SEC4 we review the Chebychev collocation method along with the Clenshaw-Curtis algorithm to compute the phase integral. In §LABEL:SEC5 we illustrate the efficiency of the combined WKB-spectral method on some numerical examples, and we conclude in §LABEL:SEC6.
2 Review of the WKB-method for the stationary Schrödinger equation
In this section we briefly review the WKB-based numerical method from  for solving the following scalar, highly oscillating initial value problem (IVP):
with possibly complex valued initial conditions. For the rest of the paper we make the following assumptions on the coefficient function :
Hypothesis A Let be a fixed smooth (real valued) function, satisfying in , which means that we are in the oscillatory regime.
Besides, let be an arbitrary real number with
This Hypothesis implies that there are positive constants , such that
and that .
Clearly, could only be piecewise . Then the ODE problem could just be restarted at an interface point of reduced regularity.
The WKB-method consists of two steps: first an analytic transformation of (LABEL:SchIVP) into a smoother problem, and then the numerical discretization of the latter.
2.1 Reformulation of the continuous problem
This transformation involves three steps: using the substitution (LABEL:U-trafo), the ODE (LABEL:SchIVP) is first transformed to a system of first order differential equations:
with the two matrices
Next the dominant part (w.r.t. ) of the resulting system matrix (LABEL:EQU), i.e. , is diagonalized by the following change of variable:
with the unitary matrix
The final transformation step eliminates the leading oscillations by using the diagonal matrix
with the -dependent, (real valued) phase function . The change of unknown
finally leads to the system
Here, the matrix
is off-diagonal, -dependent – in fact highly oscillatory, but bounded independently of .
The principal idea of the WKB-method for is as follows: instead of solving the highly oscillatory problem (LABEL:SchIVP) on a fine mesh, one solves numerically the smooth problem (LABEL:EQZ) on a coarse mesh, possibly with . Then the original solution is recovered by
As in , the above transformation assumes that the phase is known exactly on the considered interval . In special cases (like piecewise linear coefficient functions ) this is indeed possible (since and are then explicitly integrable). Then, the crucial first error terms in (LABEL:error_U), (LABEL:error_U_2ORD) would be absent. But in general, has to be obtained by a numerical quadrature, yielding an approximation . This approximate phase will be used only in the final numerical method and in our error analysis, which generalizes the analysis from .
2.2 Numerical discretization of the transformed problem
In this section we review the first () and second () order discretizations of (LABEL:EQZ). Since the matrix is highly oscillatory w.r.t. , the key issue is to find an –uniform discretization of the resulting oscillatory integrals. Let be a discretization of the interval and . In the following, for simplicity, we shall often denote the cell simply by , i.e. and .
The marching method from  is based on a truncated Picard iteration for (LABEL:EQZ) with :
with the matrices
Since is an off-diagonal matrix, the matrices are alternatingly off-diagonal (Hermitian) and diagonal (with complex conjugate entries). Their entries are (iterated) oscillatory integrals (without stationary points) involving the two small parameters and . Since the error of their numerical approximation needs to be small w.r.t. both parameters, the integrals need to be expanded w.r.t. both of them (essentially based on carefully chosen integration by parts; cf.  for a review on the numerical treatment of oscillatory integrals).
To present the discrete analogs of (LABEL:Decomp_Z_bis)-(LABEL:EXP_M) we need some notation. We define the following functions
and the phase increments
Now we recall from  the two marching methods for the vector .
First order scheme. Let be the initial condition and let . Then we define
with the matrix
Second order scheme. Let be the initial condition and let . Then we define
with the matrices
In order to compute now the numerical approximation of (LABEL:EQU) and thus to obtain the wave function as well as its derivative , we have to transform back via
In these schemes we assumed so far that the phase , and the functions and , (which involve up to five derivatives of ) are explicitly “available”. For and this is feasible, as they involve only derivatives of . But typically the phase and phase increment have to be replaced by their numerical approximates and . At this point we have two options: On the one hand we could use the approximate along with the exact functions , , and . On the other hand it appears more consistent to combine with its exact derivatives that lead to numerical approximates , , . Since we shall use a spectral integration of the phase (cf. §LABEL:SEC4), the exact derivatives of are readily available. We shall use the second option here, as this will simplify the error analysis in §LABEL:SEC3. Moreover, the error plots of these two versions are almost indistinguishable.
With the replacements , , , and , the above matrices (LABEL:1ORD_A), (LABEL:2ORD_A1), and (LABEL:2ORD_A2) shall be called , , and . Hence, the first and second order methods with approximate phase read
For both methods, the corresponding back-transformation to the variable then reads
Note that these two schemes are WKB-schemes without numerical phase approximation for the perturbed (transformed) Schrödinger equation (LABEL:EQU) with the coefficients , .
In Section LABEL:SEC3 we shall give a complete error analysis of the first order scheme (LABEL:1ORDtilde), (LABEL:Transfo_ZU_tilde) and of the second order scheme (LABEL:2ORDtilde), (LABEL:Transfo_ZU_tilde). This is a completion of the error estimates given in , since that paper used the approximate phase only in the back-transformation (LABEL:Transfo_ZU_tilde) but not in the matrices , , and .
3 Error analysis including phase errors
First we decompose the exact phase from (LABEL:PH) as with
For the approximate phase we shall make the following assumptions:
Hypothesis B Let with .
This Hypothesis implies that there are positive constants such that
with some . In the subsequent error estimates we shall use the following error bounds on that hold uniformly in :
with some positive numbers , , . In accordance with (LABEL:phi-decomp) we also define , , but we do not require that holds (cp. with (LABEL:beta-def)). We remark that this equality was also not used in the error analysis of .
Note that we assume here that is a continuous (and smooth) function on , and it is not only defined on the grid points . In particular, this is satisfied for the spectral approximation constructed in §LABEL:SEC4 below.
As a first step of the error analysis we shall estimate the error between the (continuous) solution to (LABEL:EQZ) and its perturbed analog , which is the exact solution to
with the matrix
Let the coefficient function satisfy Hypothesis A and let satisfy Hypothesis B. Then we have
with some generic constant independent of .
Proof. Step 1 (bound on the solution propagator:) To estimate the growth of the solution we compute
Using the abbreviation for , this implies for the solution propagator of (LABEL:EQZ):
and analogously for the solution propagator of (LABEL:EQZ1):
Step 2 (bound on :) We denote , which satisfies
The solution of this inhomogeneous equation reads
where we skipped in the integrands the argument ”” for brevity. In (LABEL:DZ-repr) we used the following integration by parts formula involving the propagator for some linear evolution equation :
which can be verified easily by using . Note that the integration by parts in the oscillatory integral of (LABEL:DZ-repr) will allow to recover one more -power in the estimate of . This strategy was already used in Proposition 2.2 of . In the last line of (LABEL:DZ-repr) we used also
On the r.h.s. of (LABEL:DZ-repr) we have to consider two types of differences: First we estimate
where we used for the first term in (LABEL:est1) both the trivial estimate and the mean value theorem for vector functions. We also used .
Secondly we estimate:
From (LABEL:S-est), (LABEL:tildeS-est), and (LABEL:EQZ1) we recall that the propagator , the solution , and are uniformly bounded in , , and . Hence, the estimates (LABEL:est1) and (LABEL:est2) yield the result (LABEL:Z-error) with a constant that depends only on , , , , , , , and .
This lemma allows to derive the main result of this section:
Let the coefficient function satisfy Hypothesis A and let satisfy Hypothesis B. Then the first order scheme (LABEL:1ORDtilde), (LABEL:Transfo_ZU_tilde) and the second order scheme (LABEL:2ORDtilde), (LABEL:Transfo_ZU_tilde) satisfy the following error estimates:
with some generic constant independent of , , and .
Let us compare this result with the estimates (LABEL:error_U), (LABEL:error_U_2ORD) that are due to [1, Th. 3.1]: The first error terms on the r.h.s. of (LABEL:error_Unew) and (LABEL:error_U_2ORDnew) are generalizations to -independent numerical integrations of the phase integral. The new (additional) third terms are due to using the perturbed phase in the WKB-method.
Proof. [of Theorem LABEL:ZU-est] First we estimate the error in the -variable, for :
where we used Lemma LABEL:Z-diff and the WKB-error estimate for that is analogous to (LABEL:error_U), (LABEL:error_U_2ORD) (skipping the first term, cf. [1, Th. 3.1]).
Next we estimate the error in the -variable, using (LABEL:Transfo_ZU), (LABEL:Transfo_ZU_tilde):
where we used an estimate like for the first term of (LABEL:est1), the -uniform boundedness of (see (LABEL:S-est)), and the unitarity of the matrices , . Combining the two estimates yields the result.
4 Spectral integration of the phase
The estimation of the numerical errors (LABEL:error_U) and (LABEL:error_U_2ORD) in the computation of a solution to the Schrödinger equation in the semi-classical limit via the approach of  indicates that the problematic term for small is the first on the right hand sides of these expressions. It arises from the numerical computation of the phase (LABEL:PH) and is not present if the latter can be computed exactly. In cases where this is not possible, a high order method is recommended to reduce as much as possible the role of the term proportional to . We use here spectral methods which are known to approximate analytic functions with spectral accuracy, i.e., an error decreasing exponentially with the number of modes. The numerical error for functions in accordance with Hypothesis A is known to decrease faster than any power of , which means in practice an exponential decrease, too, see e.g. . Concretely we apply a Chebychev collocation method and use the Clenshaw-Curtis  algorithm for the integral in (LABEL:PH). For points in between collocation points of the Clenshaw-Curtis algorithm, we use barycentric interpolation, see .
The basic idea of spectral methods is to approach a function on the interval via functions which are globally smooth on the considered interval. We will use here Chebychev polynomials since Chebychev series are related to Fourier series for which efficient numerical algorithms exist. Since any finite interval can be mapped via to the interval , we present all algorithms for the latter interval. We approximate via
where the Chebychev polynomials are defined as
The idea of a collocation method is to introduce collocation points , on and to impose in (LABEL:approx) equality at the collocation points,
These are equations to determine the spectral coefficients , . Choosing the as the Chebychev collocation points , , the equations (LABEL:coll) take the form
Thus the spectral coefficients are given by the discrete cosine transformation (DCT) of the function at the collocation points. Since the DCT is related to the discrete Fourier transform, it can be computed with the fast Fourier transform algorithm after some preprocessing, see for instance . Thus one advantage of a Chebychev collocation method is that a fast algorithm to compute the spectral coefficients exists.
To integrate a function approximated by the Chebychev sum (LABEL:approx), a very efficient algorithm exists due to Clenshaw and Curtis . The basis of the algorithm is the well known identity for Chebychev polynomials (simply a consequence of the addition theorems for trigonometric functions)
The antiderivative of a function approximated as a Chebychev sum (LABEL:approx) can itself be approximated by such a sum,
where the follow from the via (LABEL:chebint). In  the numerical error of the Clenshaw-Curtis algorithm was discussed showing that it is a spectral method. Identity (LABEL:chebint) can obviously also be used to approximate derivatives of functions in coefficient space. Writing , one gets from (LABEL:chebint)
i.e., the derivative of is approximated by the action of the differentiation matrix on the vector of coefficients with components , . We use these matrices to compute the derivatives appearing in the definition of the , (LABEL:def-betak). Thus these derivatives are also computed with spectral accuracy. For –error bounds of the Chebychev spectral approximation (and its derivatives) we refer to Theorem 5 and 6 in . We recall that such estimates were assumed for the error analysis in §LABEL:SEC3.
As an example we consider the function appearing in the examples in the following section for . The difference between the Clenshaw-Curtis approximation of and the exact value (the well known error function computed in Matlab to machine precision via ) in dependence of the number of collocation points is shown on the left of Fig. LABEL:ccfig. It can be seen in the semilogarithmic plot that the numerical error decreases exponentially with the number of collocation points up to where the numerical error reaches the saturation level (we work here in double precision, thus the accuracy is limited in practice to the order of because of rounding errors).
The Clenshaw-Curtis algorithm gives in principle only the antiderivative of a function at the collocation points. However in the present context, the former will be needed on more general values of . Since the basis of the approach is a sum of Chebychev polynomials, intermediate values can be obtained in principle from formula (LABEL:approx). A numerically stable and very efficient way to interpolate is to use Lagrange interpolation in the barycentric form, see  and references therein. For Chebychev collocation points, the interpolation weights can be given explicitly, and a Matlab code for this case can be found in . The difference between the antiderivative of the function and in dependence of can be seen for in Fig. LABEL:ccfig on the right. The difference on the collocation points is marked with red ‘s’. It can be seen that the error introduced by interpolation at intermediate points is also smaller than .
As mentioned above, the exponential decay of the numerical error with in the Clenshaw-Curtis algorithm for functions analytic in a strip around the real axis in the complex plane is a general feature of Chebychev series for such functions. This can be seen on the left of Fig. LABEL:ccfigN where the Chebychev coefficients of the antiderivative of are shown. It can be seen that they decrease exponentially, and that the numerical error due to truncation of the series at terms is actually due to the highest order spectral coefficient. Thus the Chebychev coefficients also provide an approach to estimate the numerical error due to a spectral method by studying the spectral coefficients. For the example , we had just considered in Fig. LABEL:ccfig the term in the integral of (LABEL:2orderWKB) since we had an independent way to compute the exact integral via the error function. This is not the case for the terms proportional to in the integral in (LABEL:2orderWKB). But the spectral coefficients of the antiderivative of these terms ( for the example ) on the right of Fig. LABEL:ccfigN indicate a similar behavior of the error as in Fig. LABEL:ccfig: for the Chebychev coefficients are of the order of the rounding error, and further increase of the number of coefficients no longer leads to higher accuracy.
The example in (LABEL:ccfig) shows that a spectral approach for functions allows in practice to reach machine precision with low resolution, here with just 14 Chebychev polynomials. Thus the numerical error reaches a plateau which itself will increase with . The latter is due to the fact that rounding errors pile up with larger values of . With finite difference methods, considerably larger values of (or equivalently smaller values of ) are needed to reach the saturation of the numerical errors, and because of this, the plateau is reached in practice at much higher values than here, of the order of (see for instance examples in ).
Note that in , the saturation level of the numerical errors was not reached since quadruple precision was used, and since the values of were not small enough to get there with the used precision. Here, we work in double precision and will reach the saturation level in most cases.
5 Numerical results
We shall present now numerical results obtained with the first and second order WKB-schemes from Section LABEL:SEC2. For our numerical tests we chose on the spatial interval with a uniform grid, and the initial condition . For both schemes we shall compare the results obtained with two versions of the numerical phase computation: on the one hand by the composite Simpson rule (with error order ) on the WKB-grid , and on the other hand by the spectral method from §LABEL:SEC4 along with barycentric interpolation at the WKB-grid points . Note that we use a Chebychev grid with points as in the previous section for the computation of the phase in (LABEL:2orderWKB) and then interpolate to the equidistant grid for the WKB-scheme. This is necessary since Chebychev collocation points are not equidistant. A consequence of this approach is that the spectral method computes the phase always to machine precision. In fact, the absolute and relative errors of is for this example always of the order .
Fig. LABEL:1ORDfig shows the results for the first order method. Plotted are the –errors of the numerical solution as a function of for 5 values of . The reference solutions were obtained with the same method, but on a much finer grid. For simplicity we shall refer in our discussion to the error estimate (LABEL:error_U). The left plot is obtained with the Simpson rule to compute . For the second term (i.e. the WKB-error) in (LABEL:error_U) dominates and the method is clearly first order in , as indicated by the upper slope triangle. also shows this behavior for . For smaller values of and large step sizes (e.g. ) the error behaves like the first error term in (LABEL:error_U), i.e. due to the Simpson rule, as visualized by the lower slope triangle. This leads to an inversion of the 5 error curves at ( at the top,