Splitting methods for the time-dependent Schrödinger equation

17 January 2000

Physics Letters A 265 Ž2000. 35–42 www.elsevier.nlrlocaterphysleta

Splitting methods for the time-dependent Schrodinger equation ¨ S. Blanes a

a,b

, P.C. Moan

b,)

Departament de Matematiques, UniÕersitat Jaume I, 12071-Castellon, ` ´ Spain b DAMTP, UniÕersity of Cambridge, Cambridge CB3 9EW, UK

Received 20 August 1999; received in revised form 11 November 1999; accepted 7 December 1999 Communicated by P.R. Holland

Abstract Cheap and easy to implement fourth-order methods for the Schrodinger equation with time-dependent Hamiltonians are ¨ introduced. The methods require evaluations of exponentials of simple unidimensional integrals, and can be considered an averaging technique, preserving many of the qualitative properties of the exact solution. q 2000 Published by Elsevier Science B.V. All rights reserved. PACS: 02.60.Cb; 02.30.Hq; 03.65.Fd Keywords: Time-dependent Schrodinger equation; Time-dependent Hamiltonian; Splitting methods; Magnus expansion; Symplectic; Lie ¨ series

1. Introduction In this work we present new fourth-order geometric integrators for solving the time-dependent Schrodinger equation Ž " s 1. ¨ i

d dt

c Ž x ,t . s Hˆ Ž p, x ,t . c Ž x ,t . ,

c Ž x ,0 . s c 0 Ž x .

Ž 1.

where c Ž x,t . is the wave function associated with the system, p ' yi EEx and Hˆ Ž p, x,t . is a Hermitian Hamiltonian operator governing the evolution of the system. Some standard techniques are found in the

)

literature for solving Ž1. as a linear ordinary differential equation: Ži. Spatial discretisation of c (x,t). Let us assume that the system is defined in the interval x g w x 0 , x f x. We can then split this interval in N parts of length ^ x s Ž x f y x 0 .rN and consider c n s c Ž x n ,t . where x n s x 0 q n^ x, n s 1, . . . , N, thus obtaining the finite dimensional linear equation

Corresponding author. E-mail: [email protected]

i

d dt

cŽ t . s HŽ t . cŽ t . ,

c Ž 0. s c 0 ,

Ž 2.

where c s Ž c1 , . . . ,c N . T g C N and H g C N=N is an Hermitian matrix associated to the Hamiltonian Žusually it is real and symmetric.. If the Hamiltonian takes the form Hˆ s TˆŽ p,t . q VˆŽ x,t ., then Vˆ is diagonal in the coordinate space and Tˆ is diagonal in momentum space. We can use complex Fast Fourier Transforms ŽFFTs. for evaluat-

0375-9601r00r$ - see front matter q 2000 Published by Elsevier Science B.V. All rights reserved. PII: S 0 3 7 5 - 9 6 0 1 Ž 9 9 . 0 0 8 6 6 - X

S. Blanes, P.C. Moan r Physics Letters A 265 (2000) 35–42

36

ing it as Tˆc Ž x n ,t . s Fy1 D T Fc Ž x n ,t ., where D T is a diagonal operator. Also finite difference methods produce similar finite dimensional systems, but in this case the spatial accuracy is reduced. Žii. Spectral decomposition. Let us suppose the solution of Ž1. can be approximated by N

c Ž x ,t . ,

Ý c n Ž t . eyi t E cn Ž x . , n

Ž 3.

preserving some qualitative properties Žsymplecticity, unitarity, time-symmetry, etc.. w3,4x. In conclusion, high-order and very efficient methods have been developed for numerically evaluating Ž4.. If the Hamiltonian is explicitly time-dependent, the situation is different. The unitary evolution oper˙ s HŽ t .U and the ator has to satisfy the equation iU exact solution can formally be written as

ns0

where En and cn are the exact eigenvalues and eigenfunctions of a time-independent Hamiltonian, Hˆ0 , e.g. a Harmonic oscillator, and the complex coefficients c n give the probability amplitude to find the system in the state cn . Then, substituting Ž3. into Ž1. we obtain a matricial equation similar to Ž2. for this basis, with c s Ž c1 , . . . ,c N . and ŽHŽ t .. i j s ² c i < Hˆ Ž t . y Hˆ0 < c j :e iŽ E iyE j .t , i, j s 1, . . . , N, where the standard bracket notation is used. In the following we will present new numerical methods of fourth-order Žin the time-step t . for solving Ž2.. They will be written in terms of exponentials of linear combinations of H0t HŽ t . dt and H0t t HŽ t . dt, preserving many of the qualitative properties of the exact solution. After this averaging of the time-dependent part of the Hamiltonian, the system can be considered as autonomous. So standard techniques like splitting methods can be used for evaluating the exponentials. To avoid the need for complex number computations, we apply the fact that the quantum system Ž2. is equivalent to a 2 N degree of freedom classical Hamiltonian system. Finally the performance of the methods is illustrated with several examples.

2. Second and fourth-order methods

t

H0

where U is the unitary evolution operator. For this problem, a number of methods have been developed. The optimal method varies for each specific problem Ždepending on the structure of H, how frequently output is desired, etc. w1x.. In most cases the geometric integrators frequently used in the context of classical mechanics w2x can be used to solve Ž4.,

H0

/,

H Ž t . ,H Ž s . ds s 0,

; t ) 0,

where wA ,Bx s AB I BA, the solution is given by

ž

U Ž t ,0 . s exp yi

t

H0 H Ž s . ds

/

.

This is not the usual case and in general one has to look for numerical approximate solutions. If we split the interval Ž0,t . into Nt subintervals then Nty1

U Ž t ,0 . s

Ł U Ž t n q t ,t n . ,

ns0

with t s trNt and t n s nt . If t is sufficiently small and HŽ t . is a bounded matrix, each step can be written in exponential form U Ž t n q t ,t n . s e V Ž t nqt ,t n ., where V s Ý`ns 1 V n is the Magnus series w5,6x. The first two terms are given by t nq t

Ht

Ž 4.

H yi t Ž s . d s

T is the standard time ordering operator which whereT gives rise to more complicated expressions to disentangle. If

V 1 s yi

If H is a constant matrix, the Schrodinger equa¨ tion has the solution c Ž t . s U Ž t ,0 . c Ž 0 . s exp Ž yitH . c Ž 0 . ,

ž U Ž t ,0 . s T e

H Ž t 1 . dt 1 ,

Ž 5.

n

V 2 s y 12

t nq t

t1

n

n

Ht Ht

H Ž t 1 . ,H Ž t 2 . dt 2 dt 1 ,

Ž 6.

and it is well known that expŽ V 1 . provides a second order approximation. Then, if the Hamiltonian is separable in the form H Ž p, x ,t . s T Ž p,t . q V Ž x ,t . ,

Ž 7.

where T and V are the matrices associated with kinetic and potential energy respectively, we can

S. Blanes, P.C. Moan r Physics Letters A 265 (2000) 35–42

evaluate Ž5. using different second order quadratures Žmidpoint, trapezoidal, etc.., e.g. a staggered grid t V 1 Ž t . s t T1r2 q Ž V0 q V1 . q O Ž t 3 . , Ž 8. 2 with T1r2 s yi TŽ t n q 2t . and Vk s yi VŽ t n q kt ., k s 0,1, and we can use second order splitting methods in several forms, e.g. U s e V 1 q O Ž t 3 . s et r2 V1 et T1r 2 et r2 V0 q O Ž t 3 . ,

Ž 9. belonging to the family of standard second order splitting methods widely analysed in the literature Žknown as leap-frog, Stormer or Verlet.. In order to ¨ see that the error is effectively of order O Žt 3 . we can utilise the BCH formula. The first error term of Ž9. is proportional to t 2 wT1r2 ,V1 y V0 x, but V1 y V0 s O Žt . if we assume sufficient smoothness of VŽ x,t .. Similar methods to Ž8. and Ž9. can be obtained by interchanging T and V. If T and V share the same time-dependent functions, it may be convenient to evaluate both using the same quadrature Žin order to evaluate the functions at the same points.. Observe that, if no output is desired, Ž9. necessitates only one time-dependent evaluation of VŽ x,t . and TŽ p,t . per step Ža FSAL property.. The matrices expŽV. are diagonal in the configuration space while expŽT. are diagonal in the momentum space. So complex FFTs can be used, requiring only O Ž N log N . floating point operations to transform between these spaces. These methods are cheap, easy to implement, have good stability properties and preserve many of the qualitative properties of the exact solution. The main problem is low accuracy. To obtain higher order methods while preserving these important properties is a more complicated task than in the case of constant Hamiltonians. Different methods have appeared for solving this problem w7–11x, but they usually are difficult to implement, are expensive or can be used only under very special conditions on the Hamiltonian. Next we will present new fourth-order numerical methods in terms of products of exponentials of linear combinations of the integrals H0t HŽ t . dt and H0t t HŽ t . dt, for a time-step t . These new methods will preserve the same qualitative properties of the exact solution as the previous second order methods.

37

It is known that expŽ V 1 q V 2 . gives a fourth order method in the time-step t w12x, but the presence of multidimensional integrals of commutators of HŽ t . can make it difficult to apply in some problems. In order to obtain methods which can be used with different quadratures, the V i of the Magnus series have been written in terms of unidimensional integrals w13x, where

V w4x ' V 1 q V 2 s H Ž0. y w H Ž0. , H Ž1. x q O Ž t 5 . ,

Ž 10 . with H Žk. s y

i

tk

t nq t

Ht

n

ž ž

s y tn q

t 2

k

//

H Ž s . ds, Ž 11 .

where k s 0,1. If we can not solve the integrals analytically Žor it is expensive., we can do it numerically. Observe that no conditions on the points where HŽ t . has to be evaluated is required. At the same time, the evaluations used for approximating H Ž0. can be reused for the calculation of H Ž1. so, we can obtain the whole method at the cost of the evaluations of only one quadrature. If HŽ t . is only known at equidistant points it is interesting to use the Simpson rule but in general the use of the Gaussian quadrature will give the most accurate scheme. In order to avoid the presence of Žusually troublesome. commutators, we will follow the technique used for splitting methods for separable Hamiltonians. Using the BCH formula and taking Ž10. for order conditions, we can obtain, for example, the following fourth-order approximations with two and three exponentials, w4x

U s eV q O Žt 5 . s e1r2 H Ž1.

Ž0.

q2 H Ž1.

Ž0.

e1r2 H Ž1.

Ž0.

y2 H Ž1.

s e H e H eyH q O Ž t 5 . ,

qOŽt 5.

Ž 12 . Ž 13 .

where we have used the fact that H Ž0. s O Žt . and H Ž1. s O Žt 2 .. It is straightforward to prove that the new methods are time-symmetric. That is: if we integrate the equations from 0 to t and next we integrate from t to 0 using the numerical method after changing the sign of the time Ž t . we will recover the same initial conditions.

S. Blanes, P.C. Moan r Physics Letters A 265 (2000) 35–42

38

A deeper analysis in order to look for more sophisticated and efficient methods for special structures of HŽ t . as well as higher order methods will be analysed in w17x. In most cases, the schemes Ž12. and Ž13. have similar accuracy. If the Hamiltonian is separable in exactly solvable parts, as is the case of Ž7. then, the linear combination a H Ž0. q b H Ž1. s Ž a T Ž0. q b T Ž1. . q Ž a V Ž0. q b V Ž1. . is also separable in two exactly solvable parts, where T Ž k . and V Ž k . have the same structure as H Ž k . after substituting Ž11. into H by T and V, respectively. Then, each exponential can be evaluated using standard splitting methods w14–16x. It is very important to notice that if in Ž7. T is quadratic in the momenta then it is possible to use Runge– Kutta–Nystrom ¨ ŽRKN. methods, which are, in general, considerably more efficient than standard Partitioned Runge–Kutta ŽPRK. methods. In w17x we will analyse the more efficient splitting methods we found in the literature in the following cases: Ž1. T and V have general structure; Ž2. T is quadratic in the momenta ŽRKN.; Ž3. both T and V are quadratic in the momenta and coordinates, respectively Žwe will refer to them as RKN2.. 3. Classical approximation In general, HŽ t . is a real symmetric time-dependent matrix and it is then possible to avoid working with complex vectors considering the following Ndimensional real vectors: q s Re c and p s Im c. The evolution of the system is then given by d q q s J m HŽ t . , p p dt with 0 HŽ t . J m HŽ t . ' . Ž 14 . yH Ž t . 0

½5

½5

ž

/

This equation is usually related to the evolution of a classical Hamiltonian system with q and p being the coordinates and momenta. The corresponding Hamiltonian function is H s 12 pT HŽ t .p q 12 qT HŽ t .q. Our approach is to apply the methods presented in the previous section to this problem. Then, for obtaining second order methods we can define 0 Vc , k ' yH k

ž

0 0 0 ; Tc , k ' 0

/

ž

Hk , 0

/

Ž 15 .

with H k s HŽ t n q kt ., k s 0, 12 ,1. Then a scheme similar to Ž9. can be used, where expŽTc, k . s I q Tc, k and expŽVc, k . s I q Vc, k . In this case it is more efficient to replace T1r2 by 12 ŽTc,0 q Tc,1 . in Ž9., hence the method will require only one evaluation of HŽ t . per step. In the calculation of the vector-matrix multiplications H q and H p we can consider real FFTs when evaluating the kinetic energy. Finally, since the last exponential can be evaluated together with the first one in the next step, the method requires one evaluation of the time-dependent matrix HŽ t . and four real FFTs per step, which is equivalent with two complex FFTs giving a cost similar to that of Ž9.. There are, however, some qualitative differences. Scheme Ž9. applied to Ž7. preserves unitarity because each exponential is an unitary operator. On the other hand in the classical case, the unitarity is preserved only if the vector Ž q, p . evolves through an orthogonal transform, which is not the case if we split the matrix in the form Ž15.. The unitarity is not preserved exactly but, in most cases it is not seriously perturbed due to the preservation of symplecticity w3,4,10x Žsymplecticity is a necessary, but not sufficient, condition for unitarity.. Scheme Ž9. is useful only if the Hamiltonian can be split in two exactly solvable parts, while scheme Ž15. can be used independently of the structure of H. Alternatively the extrapolation technique from w18x can be utilised, giving schemes with arbitrarily accurate unitarity while still retaining symplecticity. If we define now HcŽ k '

1

t

k

Jm

t nq t

Ht

n

H Ž s . ds'

ž

t

ž ž

s y tn q

0

HŽ k .

yHŽ k .

0

/

2

k

//

,

Ž 16 .

k s 0,1, we can use the fourth-order methods Ž12. and Ž13. substituting H Ž k . for HcŽ k .. In this case HcŽ k . is always separable in two exactly solvable parts, e.g. VcŽ k . '

ž

0 yHŽ k .

0 0 ; TcŽ k . ' 0 0

/

ž

HŽ k . , 0

/

Ž 17 .

S. Blanes, P.C. Moan r Physics Letters A 265 (2000) 35–42

with k s 0,1. Observing that w TcŽ0.,TcŽ1. x s w VcŽ0.,VcŽ1. x s 0 and TcŽ0. , TcŽ0. , TcŽ0. ,VcŽ0. s VcŽ0. , VcŽ0. , VcŽ0. ,TcŽ0.

s0

Ž 18 .

it can be shown that very efficient RKN2 methods can be used w3x.

4. Numerical experiments In order to appreciate the efficacy of the methods presented in this letter we will consider the Hamiltonian of a diatomic molecule with a linear time-dependent perturbation Hˆ s y

1

E2

2m E x2

q V Ž x . q xf Ž t . ,

Ž 19 .

where m is the reduced mass of the diatomic molecule. We will consider the Morse potential, 2 V Ž x . s D Ž 1 y eya x . , as a good approximation for the study of the vibrational states of a diatomic molecule w19x. D being the dissociation energy and a the length parameter. The unperturbed system has 24 bounded states with energies

ž

En s n q

1 2

/ ž

v0 y nq

1 2

2

/

v 02 4D

,

n s 0 . . . 23,

'

where v 0 s a 2 Drm . In particular we will study the HF molecule, whose parameters are: m s 1745 a.u., D s 0.2251 a.u. and a s 1.1741 a.u. For the time-dependent perturbation we will consider the following two cases: 1) Interaction with a laser field: f Ž t . s Acos Ž v t . . We will consider three cases: Ža. A s 0.0011025 a.u. and v s 0.193 v 0 Žthe Rabi frequency. according to w19x, Žb. A s 0.011025 a.u. and v s 0.9476 v 0 as in w10x, and Žc. v s 5 v 0 with A s 0.011025 a.u. 2) Collision with an atom: Av f Ž t. s , 2 cosh Ž B'v t . where v is proportional to the energy of the collision and A, B are parameters of the system. We will use A s 0.385 a.u., B s 2.5 a.u. which are consistent with the model of w20x, and we will take v s 50 v 0 .

39

As initial condition we take the ground state of the Morse oscillator, c 0 Ž x . s Rexp Ž yb x . exp Ž yg eya x . , with g s 2 Drv 0 , b s Žg y 1r2. a where R is the normalisation constant. The grid for the spatial coordinate x ranges from y0.8 to 4.32 with N s 64 with periodic boundary conditions assumed. To gauge accuracy, we consider the instantaneous mean energy of the diatomic molecule EŽ t . s ²cŽ t .< 12 P 2 q V( X )
H

Observe that now the matrices are time-independent and the problem will be to evaluate the exponentials when using the methods Ž12. and Ž13.. As mentioned previously, there myriads of methods available for evaluating the exponential of a matrix and the proper choice depends on the particular Hamiltonian. Observe that the exponentials of Ž17. are very cheap to evaluate through Žk.

eV c

Žk.

eT c

q q s ; p p y HŽ k . q

½5 ½ ½5 ½

q q q HŽ k . p s p p

5

5

Ž 21 .

requiring only one real FFTrinverse FFT for evaluating the kinetic part in H Ž k . q and H Ž k . p. Then, we can use splitting methods in the form m

et Ž AqB . s Ł e a it A e b it B q O Ž t k . ,

Ž 22 .

is1

with k ) 4. Then, for the exponentials in Ž12. we can consider that 12 HcŽ0. " 2 HcŽ1. s A q B with A s 12 VcŽ0. " 2VcŽ1. and B s 21 TcŽ0. " 2TcŽ1., where 21 H Ž0. " 2H Ž1. s 4h P 2 q 2h V( X ) q X Ž 12 f Ž0 . " 2 f Ž 1 . .. Property Ž18. allows us to use RKN2 methods, and in w3x m we can find solutions for  a i ,bi 4is1 for several values of m and k up to k s 12. These families of methods are very accurate, have good stability prop-

40

S. Blanes, P.C. Moan r Physics Letters A 265 (2000) 35–42

erties and preserve unitarity up to order k q 1. If we use one of these methods for expŽ 12 HcŽ0. q 2 HcŽ1. . and its adjoint for expŽ 12 HcŽ0. y 2 HcŽ1. . then, the total number of stages is reduced by one, leaving the total composition symmetric and in addition the unitarity is preserved up to order k q 2. In our examples we will consider the methods of w3x with m s k s 4 and m s k s 6. In case of scheme Ž13., expŽ HcŽ0. . can be evaluated in a similar fashion and expŽ HcŽ1. . is exactly and easily computed because H Ž1. is a diagonal matrix. Thus the overall cost per step will be less than for scheme Ž12.. In the following we denote schemes Ž12. and Ž13. by 2EX-k and 3EX-k respectively, with k s 4,6 depending which splitting method of w3x is applied to compute the exponential. For comparison we consider the approach presented in w10x for problem Ž14. with associated Hamiltonian H s 12 p T HŽ t .p q 12 qT HŽ t .q. This Hamiltonian can be rewritten in the form H˜ s Ž 12 p T H Ž p 2 t . p q p 1 t . q Ž 12 q T H Ž q 1 t . q y q 2 t . , where the new coordinates q1 t and q2 t and their associated momenta p 1 t and p 2 t are introduced in order to obtain a time independent Hamiltonian, which is separable in two exactly solvable parts. Observe that H is quadratic in the coordinates and momenta, but this is not necessarily the case of H˜ . Then, the number of allowed methods is considerable reduced, i.e. it is not possible to use RKN2 methods. We will consider the five stage fourth order Ž S,m s 5. and the nine stage sixth order Ž SS,m s 9. methods given in w16x, which are the most efficient fourth and sixth order symmetric PRK method we found. The methods will be referred to as McL-4 and McL-6. In Fig. 1 we consider the laser field perturbation through 100 periods of the laser frequency ŽT s 100 2vp . and display the average relative error in energy Žthe average is taken over 5–10 consecutive solution points.. The time-step used is very similar for all three examples and is chosen such that all methods require 16000, 4000 and 640 FFTs per period respectively. There are two sources of error in our methods: the truncated Magnus approximation for the time-dependent part and the splitting approximation of the exponential. It seems clear that for smooth time-dependent functions the error of Magnus is very small

Fig. 1. Average relative errors in energy through 100 periods for the laser field perturbation for: a. As 0.0011025, v s 0.193 v 0 with 16000 FFTrperiod b. As 0.011025, v s 0.9476 v 0 with 4000 FFTrperiod and c. As 0.011025, v s 5 v 0 with 640 FFTsrperiod.ŽMcL-6 is in this case unstable..

and the main error comes from the splitting method. Thus in this case one should to use the most accurate and stable splitting methods Ž2EX-6.. In Fig. 1b the system is simulated near the resonant frequency and the contribution to the total error comes equally from

S. Blanes, P.C. Moan r Physics Letters A 265 (2000) 35–42

the splitting method and the Magnus approximation, so in this case best choice is to use the cheapest Magnus scheme together with a good splitting method i.e. Ž3EX-6.. Finally, for highly oscillating Žin ’t’. functions the largest errors will originate from the Magnus expansion and thus the cheapest splitting method Ž3EX-4. will be the optimal choice ŽFig. 1c.. The Magnus approximation integrates the time dependent part more accurately than the standard splitting methods with the number of time-dependent function evaluations considerably reduced. Furthermore the accuracy improves with increase in frequency. In Table 1 we present, for the near resonant case Žas in Fig. 1b., the stability limit Žtime-step and number of FFTs., the average error in the preservation of unitarity Ž²cŽ t .
t 2EX-4 2EX-6 3EX-4 3EX-6 McL-4 McL-6

14.0 11.3 5.5 6.9 4.7 4.0

Error in FFT 700 1364 1024 1224 1500 3168

Function

unitarity y1 3

6.8=10 1.2=10y1 4 9.7=10y1 2 1.1=10y1 3 8.5=10y9 1.8=10y9

evaluations 285 182 500 333 2000 2000

41

Fig. 2. Average relative errors in energy vs number of FFTs for the collision problem, with step-sizes starting at the stability limit and from there on decreasing. Methods 2EX-4 and 3EX-4 quickly reach the machine accuracy before displaying the expected 4th order accuracy.

In conclusion we can say that the methods presented provide a superior alternative to the approach presented in Ref. w10x. The approach lends itself to high accuracy computations of Schrodinger equa¨ tions with time-dependent potentials. The methods retain symplecticity exactly while unitarity is conserved to high order Ž) 4.. In Ref. w17x we present how the approach presented in this letter for linear equations, easily can be adapted to non-linear, timedependent classical Hamiltonian systems, retaining symplecticity and time-symmetry for high accuracy, long term computation of such systems. Acknowledgements The authors thank J. Ros and F. Casas for their comments. P.C.M. would like to mention A. Iserles for suggesting the Magnus expansion for approximating oscillatory problems. S.B. acknowledges the Ministerio de Educacion ´ y Cultura ŽSpain. for a postdoctoral fellowship. The work of P.C.M. is supported by the Norwegian Research Council through contract 119089r410. References w1x C. Leoforestier, R.H. Bisseling, C. Cerjam. M.D. Feit, R. Friesner, A. Guldeberg, A. Hammerich, G. Jolicard, W.

42

w2x w3x w4x w5x w6x w7x w8x w9x

S. Blanes, P.C. Moan r Physics Letters A 265 (2000) 35–42 Karrlein, H.-D. Meyer, N. Plipkin, O. Roncero, R. Koslof, J. Comput. Phys. 94 Ž1991. 59. J.M. Sanz-Serna, M.P. Calvo, Numerical Hamiltonian Problems ŽChapman &Hall, London, 1994.. S. Gray, D.E. Manolopoulos, J. Chem. Phys. 104 Ž1996. 7099. R.I. McLachlan, S.K. Gray, Appl. Numer. Math. 25 Ž1997. 275. W. Magnus, Commun. Pure Appl. Math. 7 Ž1954. 649. P.C. Moan, J.A. Oteo, J. Ros, J Phys. A: Math. Gen. 32 Ž1999. 5133. L. Dieci, R.D. Russell, E.S. Van Vleck, SIAM J. Numer. Anal. 31 Ž1994. 261. M. Glasner, D. Yevick, B. Hermansson, Mathl. Comput. Modelling 16 Ž1992. 177. M. Hochbruck, C. Lubich, BIT. 39 Ž1999. 620.

w10x J.M. Sanz-Serna, A. Portillo, J. Chem. Phys. 104 Ž1996. 2349. w11x B.A. Shadwick, W.F. Buell, Phys. Rev. Lett. 79 Ž1997. 5189. w12x A. Iserles, S.P. Nørsett, Philosophical Trans. Royal Soc. A 357 Ž1999. 983. w13x S. Blanes, F. Casas, J. Ros, DAMTP tech. report 1999rNA08, University of Cambridge Ž1999.. w14x H. Yoshida, Phys. Lett. A 150 Ž1990. 262. w15x M. Suzuki, Phys. Lett. A 165 Ž1992. 387. w16x R.I. McLachlan, SIAM J. Sci. Comput. 16 Ž1995. 151. w17x S. Blanes, P.C. Moan, work in progress. w18x A. Iserles R. McLachlan, A. Zanna, DAMTP tech. report 1997rNA22, University of Cambridge Ž1997.. w19x R.B. Walker, K. Preston, J. Chem. Phys. 67 Ž1977. 2017. w20x B. Gazdy, D.A. Micha, J. Chem. Phys. 82 Ž1985. 4926.