26. Split-Operator Methods¶
PDF pages 1149–1166
Chapter 26 Split-Operator Methods The method of operator splitting is a general technique for making the evolution of a system numerically tractable.1 In general, the evolution (e.g., the equations of motion) of a system can be represented in terms of an operator (e.g., a matrix for a linear system). In some situations there may exist approximate decompositions of the operator into parts that are particularly simple to compute, and this is the essence of operator splitting. The unitary evolution of a quantum system according to the Schrödinger equation is particularly suited for operator-splitting methods, since the unitary time-evolution operator can often be decomposed into products of simple operators. 26.1 Splitting the Unitary Evolution Operator Consider the Schrödinger equation in one dimension for a particle in a potential V (x):
p2 2m + V (x) \psi(x, t), (26.1) where, as usual,
(26.2) We will now develop the split-operator exponential method as the basis for generating the time- dependent solutions \psi(x, t) for this Schrödinger equation. The general idea is to note that for the time- independent Hamiltonian, we can write the evolution of the wave function over a time ∆t as
(26.3) In general, the matrix exponential of H is difficult to calculate, since it is not in general diagonal in either the position or the momentum representation. However, often—as is the case for the particle Hamiltonian—the Hamiltonian is of the general form
(26.4) so that it splits into two parts, each of which is diagonal in either the position or the momentum basis. From the Baker–Campbell–Hausdorff expansion2 for arbitrary operators A and B, eAeB = exp
2[A, B] + 1 12[A, [A, B]] + 1 12[[A, B], B] + . . . , (26.5) we can see that
(26.6) 1Hans De Raedt, ‘‘Product Formula Algorithms for Solving the Time Dependent Schrödinger Equation,’’ Computer Physics Reports 7, 1 (1987) (doi: 10.1016/0167-7977(87)90002-5). 2R. M. Wilcox, ‘‘Exponential Operators and Parameter Differentiation in Quantum Physics,’’ Journal of Mathematical Physics 8, 962 (1967) (doi: 10.1063/1.1705306).
Chapter 26. Split-Operator Methods In fact, a symmetrized splitting is even more accurate:
(26.7) Note that the errors here are regarded as asymptotic, since there is no guarantee that the BCH expansion converges. These factorizations (the second, more accurate one being our preference) imply the following numerical method to evolve the wave function \psi(x, t) to \psi(x, t + ∆t):3 1. Compute the half-step for the spatial part:
(26.8) The exponential operator is diagonal in the position representation, so this is easy. 2. Compute the Fourier transform \phi(p) of the new wave function:
(26.9) In view of fast-Fourier-transform (FFT) numerical algorithms, this operation can be done efficiently. 3. Now apply the ‘‘drift’’ part of the evolution operator:
(26.10) The exponential operator is diagonal in the momentum representation, so this is also easy. 4. Compute the inverse Fourier transform \psi(x) of the updated wave function:
(26.11) 5. Finally, compute the last half-step for the spatial part, now that we are back in the representation where it is diagonal:
(26.12) This procedure is then repeated for many time steps ∆t until the desired final time is reached. The local error in each step is O(∆t3), but because the number of steps until the final time is O(∆t−1), the global error of the final solution is O(∆t2). Note that this type of symmetrized splitting works for more general splittings. For example, the decomposition
(26.13) leads to the operator splitting
(26.14) It can be shown that this type of symmetric operator splitting has O(∆t3) error for decompositions of H into any finite set of terms Hj, and the overall accuracy is still preserved even if we can only approximately compute each Hj with an error of O(∆t3). 3This type of split-operator Fourier transform algorithm was first proposed by J. A. Fleck, Jr., J. R. Morris, and M. D. Feit, ‘‘Time-Dependent Propagation of High Energy Laser Beams through the Atmosphere,’’ Applied Physics 10, 129 (1976); M. D. Feit, J. A. Fleck, Jr., and A. Steiger, ‘‘Solution of the Schrödinger Equation by a Spectral Method,’’ Journal of Computational Physics 47, 412 (1982) (doi: 10.1016/0021-9991(82)90091-2).
26.2 Time-Dependent Potentials 26.2 Time-Dependent Potentials To handle a time-dependent potential V (x, t), the time-evolution operator has the same form except for the replacement V (x)∆t −\rightarrow Z t+∆t t dt′ V (x, t′), (26.15) since we recall from Section 4.2.1 that the unitary time-evolution operator (from time t to t + ∆t) now has the form
" −i ¯h Z t+∆t t dt′ H(t′) # , (26.16) if we assume the Hamiltonian commutes with itself at different times. The simplest approximation is to take Z t+∆t t
(26.17) where we may take \tau to be any time in the interval [t, t+∆t] without changing the order of the approximation
then the factorized unitary evolution operator becomes exp " −i ¯h Z t+∆t t dt′ H(t′) #
(26.18) It is better to use an approximation for the integral of V (x, t) accurate to second order, as in the trapezoidal rule: Z t+∆t t
2 + O(∆t3). (26.19) This expression then gives the more accurate splitting exp " −i ¯h Z t+∆t t dt′ H(t′) #
(26.20) However, it turns out to have the same order of accuracy if we instead use the simpler splitting exp " −i ¯h Z t+∆t t dt′ H(t′) #
(26.21) where we evaluate the potential at the beginning and end points of the time interval in the two spatial half-steps. This splitting follows from the BCH expansion, which implies that ehA1ehBehA2 = exp
2 [A1 −A2, B] + h2 2 [A1, A2] + O(h3)
(26.22)
\tau = ∆t/2 in the above first-order method above corresponds to the midpoint rule of integration: Z t+∆t t
(26.23) Thus, we may also take exp " −i ¯h Z t+∆t t dt′ H(t′) #
(26.24) or exp " −i ¯h Z t+∆t t dt′ H(t′) #
(26.25) as symmetric, second-order splittings.
Chapter 26. Split-Operator Methods 26.3 Richardson Extrapolation Now we will introduce the generally applicable method of Richardson extrapolation4 to develop yet higher-order approximations. First, notice again that the global error associated with the splitting (26.21) is O(∆t2), one order lower than that of the local truncation error. We then take Richardson’s ansatz that the global error of the numerical solution takes the form
(26.26) where \psi(x, t) is the exact solution, ˜\psi∆t(x, t) is the numerical approximation generated with time steps of ∆t, and the en(x, t) are functions independent of the time step ∆t. Note that only the even-order terms appear here; we leave it as an exercise to give an argument for this, but it is important to note that this expansion is valid for evolution via the second-order, symmetric splitting above. Thus, by computing the solution at time t with multiple step sizes (say, ∆t, ∆t/2, ∆t/3), we can devise the proper linear combinations of these solutions that will cancel the error terms, and thus generate methods of higher global order. For example, a fourth-order method follows from using step sizes of ∆t and ∆t/2,
(26.27) The error term e2 (which has been cancelled in the method here) can also be estimated by h
i
(26.28) A sixth-order method follows from using step sizes of ∆t, ∆t/2, and ∆t/3: ˜\psi∆t(x, t) −16
(26.29) in which case the smallest cancelled error term e4 is estimated by −13
(26.30) These methods follow from writing out the Richardson expansions for \psi(t + ∆t) in powers of ∆t (for all the step sizes), setting up the equations as a linear system, and then solving them for the ‘‘true’’ solution and the error coefficients. Higher order methods are, of course, possible, but at some point become too computationally expensive to be worthwhile. Thus, we can construct a fourth-order [i.e., global error of O(∆t4)] method as follows: 1. Start with the simulated state ˜\psi(x, t) at time t. 2. Compute ˜\psi(x, t + ∆t) using a single step of ∆t, using the O(∆t3) operator splitting (26.21). Call this
- Compute ˜\psi(x, t+∆t) using a two steps of ∆t/2, using the same operator splitting. Call this ˜\psi∆t/2(x, t+ ∆t). 4. Combine the two results as in Eq. (26.27) to obtain the updated state:
(26.31) Then iterate this method until the desired final time. The global error of the solution will then be O(∆t4). The sixth-order method is a straightforward generalization of this procedure (the same is true of higher-order methods). The advantage of the higher-order methods, of course, is the possibility for improved accuracy for a given step size, and thus a larger step size (and, ideally, less work) for a given accuracy goal. On the other hand, the second-order method is unitary (i.e., it explicitly preserves the norm of \psi), while the fourth- and sixth-order methods are not. 4L. F. Richardson, ‘‘On the Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving Differ- ential Equations, with an Application to the Stresses in a Masonry Dam,’’ Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 83, 335 (1910).
26.3.1 Numerical Test¶
26.3 Richardson Extrapolation 26.3.1 Numerical Test To illustrate these methods, we will try them out on the harmonic oscillator, with Hamiltonian H = p2 2 + x2 2 , (26.32) using a Gaussian wave packet with initially \sigmax = \sigmap = 1/ \sqrt
wave packet is initially x0 = 0 and p0 = 20. The grid consisted of 1024 points, with a range in momentum of \pmpmax = \pm16 \sqrt 2\pi and a grid spacing ∆p = \sqrt 2\pi/32 (with the same range and spacing for the position
wave packet, and according to the exact solution these should be exactly equal. Below is plotted the rms error as a function of the step size ∆t, defined as
N v u u t N X j=1
, (26.33) where \phi(p) is the exact momentum-space wave function, \phij ≡\phi(pj) is the exact result on the grid, where pj is the momentum of the jth grid point, ˆ\phij is the numerical approximation to \phij, and N is the number of grid points. -5.4 Dt -1 -2 -3 -4 -5 -14 rms error -12 -10 -8 -6 -4 -2 2nd order 4th order 6th order The results for methods of orders 2, 4, and 6 are plotted along with lines to indicate scaling according to ∆t2, ∆t4, and ∆t6. From the plot, it is easy to see that the second-order method scales nicely as ∆t2 until the error is quite large. The fourth-order method also scales as ∆t4 once the error is sufficiently small. However, there is a strong deviation for large step sizes, where it is no longer possible to ignore the higher-order terms in the Taylor expansions, and the Richardson steps are far from unitary. For sufficiently small steps, the fourth-order scaling also breaks down as rounding errors dominate the rms error (note that the rounding error increases as the number of steps increases). For the same step, the fourth-order method also shows far better accuracy than the second-order method. The sixth-order method is qualitatively similar to the fourth-order method, with the same instability and roundoff behaviors, but there seems to be no regime in which the scaling is ∆t6. Strangely, the scaling seems to be somewhat faster. Note that the sixth-order method is correspondingly more accurate in the regime of ‘‘sensible’’ step sizes, and achieves a much better accuracy at the ‘‘sweet spot’’ of best accuracy compared to fourth-order. Thus we see that for a given step size, a higher-order method performs much better than a lower-order counterpart, provided the step size is small enough. However, a higher-order step costs more than a lower- order step in terms of computer time, simply because it involves more substeps. Plotted below is the same performance data, but now the error is plotted as a function of the execution time on a computer.
26.4.1 Hierarchical Construction¶
Chapter 26. Split-Operator Methods -1.5 cpu time -1 -14 rms error -12 -10 -8 -6 -4 -2 2nd order 4th order 6th order Obviously, the execution time will vary widely depending on problem-specific optimizations, language and compilers, and other details. For this test, the methods were implemented in the Octave language, with no special effort to optimize the methods (i.e., substeps that could have been combined into larger steps to save time were not combined, and so on), to give the highest-order methods the least ‘‘advantage.’’ The cpu times shown are nominally in seconds, but again this would vary widely depending on the hardware, language, and compiler used. However, it is clear that for a given (high) accuracy, the savings in cpu time in using a high-order method are enormous, potentially measured in orders of magnitude. 26.4 Unitary Evolvers It is possible to construct higher-order, explicitly unitary split-operator methods. Let’s start with the case of a time-independent Hamiltonian. Then we can write the second-order approximation for the unitary evolution operator U(t, t + ∆t) ≡U(∆t) as
(26.34) Then the idea is to look for ordered products that give order m approximations to the unitary evolution operator:
Y j
(26.35) That is, to search for products that give better approximations (m > 2) to the evolution operator. Clearly, these approximations are still unitary. 26.4.1 Hierarchical Construction There are multiple strategies for doing this, but one particularly useful strategy is the following hierarchical construction.5 The construction is that given the (m −1)th approximant ˜U (m−1)(∆t), we can construct the mth approximant by choosing
r Y j=1 ˜U (m−1)(sj ∆t), (26.36) 5Masuo Suzuki, ‘‘Fractal Decomposition of Exponential Operators with Applications to Many-Body Theories and Monte Carlo Simulations,’’ Physics Letters A 146, 319 (1990) (doi: 10.1016/0375-9601(90)90962-N); G. Dattoli, L. Giannessi, M. Quattromini, and A. Torre, ‘‘Symmetric decomposition of exponential operators and evolution problems,’’ Physica D 111, 129
26.4 Unitary Evolvers provided the step factors sj satisfy the two conditions r X j=1 sj = 1, r X j=1 s m j = 0. (26.37) The first condition assures that the total effect of the approximant is to evolve the system by a time step ∆t, since the approximant is the composition of r total substeps. The second condition ensures the cancellation of the order m error terms, which we can see by noting the form
r Y j=1 exp −iH ¯h sj ∆t
(26.38) of the evolution operator, and then substituting the (m −1)th approximant exp −iH ¯h sj ∆t
(26.39) where gm is an undetermined function representing the uncontrolled error at the mth order. We thus conclude that
r Y j=1 ˜U (m−1)(sj ∆t) + gm(H/¯h) r X j=1
(26.40) in which case the order m error term vanishes if we require the conditions (26.37) to be satisfied. In fact, recalling from the above argument that the global error expansion involves only even powers of ∆t, it similarly follows that if U(∆t) and the ˜U (m)(∆t) are unitary, and the splittings are time-reversal symmetric,
(26.41) then the local truncation error of ˜U (m)(∆t) involves only odd powers of ∆t. Thus, the recursive method steps two orders at a time, so long as the decompositions are symmetric:
r Y j=1 ˜U (m−1)(sj ∆t). (26.42) This is true providing that (26.37) are satisfied as well as the symmetry condition
(26.43) for all j. To construct a specific method we can use the splitting
(26.44) where s satisfies
(26.45) Requiring s to be a real root gives s = 2 − 2m−1\sqrt 2. (26.46) Thus, a fourth-order method is6
(26.47) 6Michael Creutz and Andreas Gocksch, ‘‘Higher-order hybrid Monte Carlo algorithms,’’ Physical Review Letters 63, 9 (1989) (doi: 10.1103/PhysRevLett.63.9). Haruo Yoshida, ‘‘Construction of higher order symplectic integrators,’’ Physics Letters A 150, 262 (1990) (doi: 10.1016/0375-9601(90)90092-3); Etienne Forest and Ronald D. Ruth, ‘‘Fourth-Order Symplectic Integra- tion,’’ Physica D 43, 105 (1990); J. Candy and W. Rozmus, ‘‘A symplectic integration algorithm for separable Hamiltonian functions,’’ Journal of Computational Physics 92, 230 (1991) (doi: 10.1016/0021-9991(91)90299-Z); André Bandrauk and Hai Shen, ‘‘Improved exponential split operator method for solving the time-dependent Schrödinger equation,’’ Chemical Physics Letters 176, 428 (1991) (doi: 10.1016/0009-2614(91)90232-X).
Chapter 26. Split-Operator Methods
s4 = 2 − 3\sqrt 2 \approx 1.3512 0719 1959 6576 3405 1 −2s4 = − 3\sqrt 2 − 3\sqrt 2 \approx −1.7024 1438 3919 3152 6810. (26.48) Similarly, iterating this procedure produces the sixth-order method
(26.49)
s6 = 2 − 5\sqrt 2 \approx 1.1746 7175 8089 3633 8450 1 −2s6 = − 5\sqrt 2 − 5\sqrt 2 \approx −1.3493 4351 6178 7267 6899. (26.50) Of course, it is straightforward to compute the coefficients for the higher-order methods. Let’s compare these methods to the Richardson-extrapolation methods. The fourth-order Richardson method requires three applications of ˜U (2)(∆t) for various step sizes for either the Richardson or unitary method. For the particular decomposition (26.4) of the Hamiltonian, note that some of the potential op- erators can be combined in both cases, but the number of Fourier transforms is the same. The sixth-order Richardson method, on the other hand, requires six applications of ˜U (2)(∆t) for the Richardson method, compared to nine applications for the unitary method. Thus, at higher orders, unitarity comes at the expense of extra computational effort. Unitarity may well be a desireable property of an integration method, since it preserves a known invariant. However, one may take the opposite view, and advocate a method that does not preserve unitarity, since the computed norm of the wave function can act as a diagnostic for the overall accuracy of the solution. On the other hand, at high orders, the methods may become quite unstable for large step sizes, since the assumption that the error terms are small breaks down. 26.4.1.1 High-Order Methods with Reduced Substep Intervals One potential concern of the methods based on the unitary three-term recursion (26.44) is that the step-size factors s2m and 1 −2s2m, by examination of Eq. (26.46), can be seen to always be ‘‘oversteps,’’ s2m > 1 1 −2s2m < −1, (26.51) so that the middle ‘‘backwards’’ step is necessary to compensate for the forward steps. The cancellation required by these large steps may cause problems with stability and accuracy of the solutions. It is possible to find other hierarchical constructions that have smaller steps, but at the expense of additional computational effort to achieve the same order. For example, the five-factor symmetric decomposition7
(26.52) is of the same form as Eq. (26.40), and thus s satisfies
(26.53) or explicitly, s = 4 − 2m−1\sqrt 4. (26.54) 7G. Dattoli, op. cit.
26.4 Unitary Evolvers so that the step factors s and 1 −4s satisfy the bounds 3 < s2m < 1 |1 −4s2m| < 2 3, (26.55) Specifically, the fourth-order values are s4 = 4 − 3\sqrt 4 \approx 0.4144 9077 1794 3757 3714 1 −4s4 = − 3\sqrt 4 − 3\sqrt 4 \approx −0.6579 6308 7177 5029 4856, (26.56) and the sixth-order values are s6 = 4 − 5\sqrt 4 \approx 0.3730 6582 7733 2728 2478 1 −4s6 = − 5\sqrt 4 − 5\sqrt 4 \approx −0.4922 6331 0933 0912 9910. (26.57) Again, at fourth order, the method here requires five applications of ˜U (2)(∆t), compared to three for either the Richardson or the three-factor unitary method. At sixth-order, the method here requires 25 unitary factors, compared to nine for the three-factor construction and six for the Richardson method. Going to higher orders by this route is computationally rather expensive, so the benefits must be carefully weighed against the extra cost in terms of time.8 26.4.1.2 High-Order Minimal-Product Methods In defining the sixth-order unitary methods above, we saw that their computational effort went far beyond that of the Richardson method. It is also possible to skip any hierarchical construction and search for high- order methods directly in terms of products of ˜U (2)(∆t), and thereby generate methods that use products of fewer operators. This is possible because there are fewer constraints on the relative time steps within each factor. This is done by directly computing the error terms in the BCH expansion. For example, seven-factor, sixth-order symmetric splittings of the form
can be constructed, but in general the equations to be solved to determine the wj factors are complicated, and the solutions must be found numerically. One such method has9 w1 \approx −1.1776 7998 4178 87 w2 \approx 0.2355 7321 3359 357 w3 \approx 0.7845 1361 0477 560, (26.59)
extrapolation of the same order; evidently, the preservation of unitarity still demands some additional work. Additionally, 15-factor, eighth-order symmetric splitting of the form
(26.60) 8There may also be advantages to further requiring that the steps all be positive. See Siu A. Chin and C. R. Chen, ‘‘Fourth order gradient symplectic integrator methods for solving the time-dependent Schödinger equation,’’ Journal of Chemical Physics 114, 7338 (2001) (doi: 10.1063/1.1362288). 9Haruo Yoshida, op. cit.; Yoshida gives two other sixth-order methods in this reference.
Chapter 26. Split-Operator Methods can be constructed, but in general the equations to be solved to determine the wj factors are complicated, and the solutions must be found numerically. An example is10 w1 \approx 0.3117 9081 2418 427 w2 \approx −1.5594 6803 8214 47 w3 \approx −1.6789 6928 2596 40 w4 \approx 1.6633 5809 9633 15 w5 \approx −1.0645 8714 7891 83 w6 \approx 1.3693 4946 4168 71 w7 \approx 0.6290 3065 0210 433, (26.61)
to the 10 required by Richardson extrapolation of the same order. 26.4.1.3 High-Order Treatment of Time-Dependent Potentials The same strategies above for unitary, high-order product methods apply to cases of time-dependent poten- tials, so long as we again make the generalization V (x)∆t −\rightarrow Z t+∆t t dt′ V (x, t′), (26.62) in each factor ˜U (2)(∆t) −\rightarrow ˜U (2)(t + ∆t, t) (26.63) that composes the higher order methods, so long as we define the time interval for each factor appropriately.11 In fact, we can use the splitting (26.21) as our fundamental second-order factor to handle the time-dependent potentials, since it has the right order of accuracy and possesses time-reversal symmetry:
(26.64) Thus, for example, the fourth-order symplectic method (26.47) becomes
(26.65) so that the fourth-order operator is the composition of the three second-order operators, and the time intervals covered by each second-order operator is the same as before. The other methods above generalize in the same way. 26.4.1.4 Numerical Test: Fourth Order Using the same test problem as for the Richardson-extrapolation test, we can compare different fourth-order methods above: the Richardson-extrapolation method, the unitary hierarchical method of Section 26.4.1, and the unitary reduced-step method of Section 26.4.1.1. 10Haruo Yoshida, op. cit.; Yoshida gives four other eighth-order methods in this reference. 11André Bandrauk and Hai Shen, ‘‘Exponential split operator methods for solving coupled time-dependent Schrödinger equations,’’ Journal of Chemical Physics 99, 1185 (1993) (doi: 10.1063/1.465362).
26.4 Unitary Evolvers -4 Dt -1 -2 -3 -4 -13 rms error -12 -10 -8 -6 -4 -2 4th order reduced step Richardson hierarchical As a function of step size, the Richardson and reduced-step methods are clearly superior to the hierarchical method in this test, with the reduced-step unitary method having the highest accuracy. Of course, both unitary methods do not suffer from the instability for large step sizes, but this is in the regime of poor accuracy anyway. The Richardson method is the best of the methods in terms of roundoff error. However, there is a wide range of complexity among the different methods, so it is useful to consider the accuracy as a function of cpu time. -1 3.3 cpu time -13 rms error -12 -10 -8 -6 -4 -2 4th order reduced step Richardson hierarchical Perhaps counterintuitively, despite the extra complexity, the reduced-step method does quite well, with the Richardson method narrowly being the most efficient in terms of cpu time for high accuracy. 26.4.1.5 Numerical Test: Sixth Order Similarly, we can compare the various sixth-order algorithms: Richardson-extrapolation, the unitary hierar- chical method of Section 26.4.1, the unitary reduced-step-size method of Section 26.4.1.1, and the unitary minimal-product method of Section 26.4.1.2.
26.4.2 Nonlinear Schrödinger Equations¶
Chapter 26. Split-Operator Methods -3 Dt -1 -2 -14 rms error -12 -10 -8 -6 -4 -2 6th order reduced step Richardson hierarchical minimal steps The Richardson method is notable here as, oddly, the only one that does not scale at sixth order, but also achieves the best overall accuracy at the optimal step size. The reduced-step method clearly achieves the largest step size for a given accuracy of all the methods. However, recall that with 25 second-order substeps, it is also by far the most computationally expensive on each step, so it is important to consider cpu time for any target accuracy. -0.7 cpu time -14 rms error -12 -10 -8 -6 -4 -2 6th order reduced step Richardson hierarchical minimal steps However, we see that the reduced-step method compares quite favorably with the other unitary methods. Despite its complexity, it seems to be a reasonable choice. However, for high accuracies it still seems to be edged out by the Richardson method. It could be that in more complicated problems, the stability properties of the reduced-step unitary method are beneficial enough to make it the method of choice, although Richardson extrapolation certainly seems to be a generally useful and simple method. 26.4.2 Nonlinear Schrödinger Equations One can also consider nonlinear Schrödinger equations with potentials of the form
(26.66) which lead to cubically nonlinear terms in the Schrödinger equation. The Gross–Pitaevskii equation is precisely of this form, and models the evolution of a Bose–Einstein condensate: the nonlinear term models
26.4.3 Symplectic Integration¶
26.4 Unitary Evolvers energy shifts due to two-body collisions, so that there is an extra spatial potential proportional to the condensate number density. Clearly, we can adapt the above methods if we take the order-two operator to be
The problem here is that we need to know \psi(x, t+∆t) to construct the leftmost factor, but that’s what we’re trying to calculate. Thus, we have an implicit method. The simplest way to handle this is to use functional iteration (or perturbation theory, if you want to think of it that way): simply compute \psi(x, t + ∆t), but using as a guess the partially evolved wave function
(26.68) Then use the result as a guess to recompute \psi(x, t + ∆t), and so on until the procedure converges. This strategy is certainly ‘‘safe’’ in adapting the above procedures, since the nonlinear potential amounts to an ordinary time-dependent potential, so long as we can calculate the potential at any time. However, things turn out to be a bit easier than this, at least at low orders.12 Observe that
(26.69) Since in the methods above, |\psi(x, t)|2 is always multiplied by ∆t2, so with an error of O(∆t3), we can evaluate \psi(x, t) in the nonlinear potential at any other time up to ∆t away. What this means, is that even with a single iteration—that is, using \psiguess(x, t) in place of \psi(x, t+∆t) in (26.67)—the order of the method is not changed. (Recall that this was not true of the ordinary potential V (x, t), which had to be evaluated at the proper points.) Thus, we have an explicit, second-order method for evolving the GPE. However, the explicit method here seems to lack time-reversal symmetry, and so it is not obvious how useful it is to compose higher-order methods. An intermediate strategy between the explicit and implicit methods is to iterate the implicit scheme once, for a Runge–Kutta style approximation to a symmetric splitting. (Note that iterating an implicit method a fixed number of times is in fact an explicit scheme, but may approximate the implicit scheme well.) However, it has been shown13 that the simple explicit scheme above works to construct a fourth-order method. It is not clear if the simple explicit strategy works at higher orders, and if higher-order methods are to be used, it is worth checking the sensitivity of the solution to the number of iterations. 26.4.3 Symplectic Integration The high-order unitary methods above are also important in classical mechanics, when simulating trajectories of Hamiltonian systems. Rather than conserving the norm of a wave function, the methods applied in this context conserve phase-space volumes, and in the case of autonomous Hamiltonian systems, they do a better job of conserving total energy.14 For example, such methods are important in celestial mechanics, where it is desireable to obtain highly accurate, long-term solutions to Newton’s equations without spurious damping effects. Recall from Section 4.4.3.2 that Hamilton’s equations can be written
\partial H
, (26.70) where {f, g}P is the Poisson bracket, the 2N canonically conjugate coordinates are combined as
(26.71) 12A. D. Bandrauk and Hai Shen, ‘‘High-order split-step exponential methods for solving coupled nonlinear Schrödinger equations,’’ Journal of Physics A: Mathematical and General 27, 7147 (1994) (doi: 10.1088/0305-4470/27/21/030); Juha Javanainen and Janne Ruostekoski, ‘‘Symbolic calculation in development of algorithms: split-step methods for the Gross- Pitaevskii equation,’’ Journal of Physics A: Mathematical and General 39, L179 (2006) (doi: 10.1088/0305-4470/39/12/L02). 13Juha Javanainen and Janne Ruostekoski, op. cit. 14For further reading and the source of some of the notation here, see Denis Donnelly and Edwin Rogers, ‘‘Symplectic integrators: An introduction,’’ American Journal of Physics 73, 938 (2005) (doi: 10.1119/1.2034523).
Chapter 26. Split-Operator Methods and Ω\alpha\beta is the canonical cosymplectic form,
0n In −In 0n , (26.72) with In denoting the n \times n identity matrix and 0n the n \times n null matrix. Suppose we write the Poisson bracket as a differential operator,
(26.73) where for any phase-space function f(z),
\partial H
. (26.74) Noting that DH acts individually on each component of z\alpha, we can simply write the vector relation
(26.75) The formal solution of Eq. (26.75) is in terms of the exponentiated operator:
(26.76) Now we will specialize to one degree of freedom, where the generalization to higher dimensions is simple. If the Hamiltonian splits into the form
(26.77) as before, then the Poisson-bracket operator takes the explicit form
(26.78)
general. We can thus use exactly the above theory to split the exponential operator exp(DH ∆t) for evolution over a short time ∆t into simpler parts of the form exp(DT ∆t) and exp(DV ∆t), where the effects of each component operator can be written explicitly. Thus, the partial kinetic solution is
x(t) p(t) = x(t) + T ′[p(t)] ∆t p(t) , (26.79) while the partial potential solution is
x(t) p(t) = x(t) p(t) −V ′[x(t)] ∆t . (26.80) Since we can write out the explicit effect of each operator on z(t), there is no need to implement Fourier transforms between operators, as in the quantum-mechanical case. However, the splittings above still provide numerical methods of the same orders as before. 26.4.3.1 Euler–Cromer Method The simplest splitting we can use is the first-order splitting (26.6), which here bcomes
(26.81) Writing this out explicitly,
x(t) p(t)
x(t) + T ′[p(t + ∆t)] ∆t p(t) −V ′[x(t)] ∆t + O(∆t2) (26.82)
26.4 Unitary Evolvers Written out as a mapping, this method becomes¶
(26.83)
see that this method is first-order accurate and half-implicit, but it easy to iterate since you just compute the new pj+1 before you compute the new xj+1. This mapping is a canonical transformation, since the transformation zj −\rightarrow zj+1 preserves the cosymplectic form:
\partial (zj)µ Ωµ\nu
. (26.84) This follows from writing out the Jacobian matrix
\partial (zj) = 1 −T ′′(pj+1)V ′′(xj) ∆t2 T ′′(pj+1) ∆t −V ′′(xj) ∆t (26.85) and multiplying out the matrix product to see that Ω\alpha\beta is unchanged: 1 −T ′′(pj+1)V ′′(xj) ∆t2 T ′′(pj+1) ∆t −V ′′(xj) ∆t −1 1 −T ′′(pj+1)V ′′(xj) ∆t2 −V ′′(xj) ∆t T ′′(pj+1) ∆t = −1 . (26.86) Also, directly from the fact that the Jacobian determinant is unity,
\partial (zj) = det 1 −T ′′(pj+1)V ′′(xj) ∆t2 T ′′(pj+1) ∆t −V ′′(xj) ∆t = 1, (26.87) we see explicitly the preservation of phase-space areas. Since these maps correspond to canonical trans- formations, they are symplectic map and thus this approximation scheme is a symplectic integrator. Compared to the Euler method, for approximately harmonic oscillation, the Euler–Cromer method produces bounded energy errors, while the Euler method produces asymptotically growing errors. We can obviously interchange the order of the operators in (26.82) to obtain another first-order sym- plectic algorithm. The mapping for this method analogous to Eqs. (26.83) is
(26.88) and it follows in the same way as the Euler–Cromer method that this mapping is symplectic. 26.4.3.2 Verlet Method It is thus straightforward to construct a second-order symplectic method by using the symmetric splitting (26.7), which adapted to symplectic integration becomes
(26.89) The mapping for this splitting can be rewritten as eDH∆t = h eDV ∆t/2eDT ∆t/2i h eDT ∆t/2eDV ∆t/2i + O(∆t3), (26.90) and thus corresponds to a composition of an Euler–Cromer step as in Eqs. (26.83) with a ‘‘reverse’’ Euler– Cromer step as in Eqs. (26.88), both of step size ∆t/2. The composition of symplectic maps is still symplectic, and thus we have a second-order symplectic method.
Chapter 26. Split-Operator Methods To write out the mapping explicitly, first we write out the composition of the two maps (26.83) and (26.88), each for a half step:
(26.91) Eliminating the intermediate step, we find
∆t. (26.92)
∆t. (26.93) This is the velocity form of the Verlet method, or the velocity Verlet method. The original form for the Verlet method15 follows from using the other possible form of the symmetric splitting, where the ‘‘drift’’ operator is applied first and last. This amounts to the replacements x −\rightarrow p, p −\rightarrow x, T −\rightarrow −V , and V −\rightarrow −T in Eqs. (26.92):
∆t
(26.94)
∆t
(26.95) Notice that in both cases, each variable is evolved according to a method symmetric about the middle time ∆t/2, either by computing the average function value (fj + fj+1)/2, which approximates the middle value fj+1/2, or by evaluating a function after taking half of an Euler-type step. The second-order nature of these approximations follow from the two second-order integral approximations Z t+∆t t
∆t + O(∆t3), (26.96) which both appear in the Verlet methods. 26.4.3.3 Higher Order Methods The nice thing about symplectic methods here is that we have already worked out the theory to extend the second-order methods to higher order. Because of the time-reversal symmetry noted above for the Verlet method, the local truncation error only involves odd powers of ∆t (or the global error involves only even 15Loup Verlet, ‘‘Computer ‘Experiments’ on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules,’’ Physical Review 159, 98 (1967) (doi: 10.1103/PhysRev.159.98).
26.4 Unitary Evolvers powers of ∆t). All the techniques above for the split-operator Fourier-transform method thus work to extend the Verlet method, if we replace each application of ˜U (2)(s∆t) by a Verlet step (in either form, as long as we stick to the same form) of s∆t. It is not useful to write out the explicit algorithms here, since in general one would simply code a subroutine to take a Verlet step and then call it multiple times on each time step of ∆t—say, three times for the fourth-order method with step sizes s4∆t, (1 −2s4)∆t, and s4∆t corresponding to the three-step method of Eqs. (26.47). 26.4.3.4 Time-Dependent Potentials In principle, the above assumption of a time-independent Hamiltonian is not restrictive, since any Hamilto- nian system of N degrees of freedom with explicit time dependence is formally equivalent to a Hamiltonian system of N+1 degrees of freedom with a time-independent Hamiltonian. (Thus, an explicitly time-dependent system of N degrees of freedom is often said to have ‘‘N 1 2 degrees of freedom.’’) However, it is useful to generalize the above methods for time-dependent Hamiltonians anyway as a more direct method for handling explicit time dependence. To handle explicitly the case of time-dependent potentials, we make the same ad- justments as for the unitary methods to evaluate the potential at the proper time. Thus, the Euler–Cromer method (26.83) is not sensitive to the precise time of evaluation, and for example becomes
(26.97) For the Verlet methods, we can take in analogy with Eq. (26.21) the splitting eDH∆t = h eDV (t+∆t)∆t/2eDT ∆t/2i h eDT ∆t/2eDV (t)∆t/2i + O(∆t3), (26.98) where DV (t) has the potential evaluated at time t, so that the velocity Verlet method (26.92) becomes
∆t, (26.99) and the usual Verlet method becomes
∆t
(26.100) preserving the second-order accuracy of the methods in the explicitly time-dependent case.
Chapter 26. Split-Operator Methods 26.5 Exercises Problem 26.1 Verify the factorizations of Eqs. (26.6) and (26.7) using the BCH formula. Problem 26.2 Justify the following statement: despite the different local errors, the error in the evolution of the quantum state according to the two splittings, Eqs. (26.6) and (26.7), after many steps is approximately the same. Problem 26.3 Why do only the even-order error terms appear in the Richardson expansion (26.26)? Problem 26.4 Verify that the Richardson-extrapolation formulae (26.27), (26.28), (26.29), and (26.30) are correct and have the advertised accuracy.