1	Higher Order Symmetric Variable Stepsize Methods. Benedict Leimkuhler, Christopher Sweet. Department of Mathematics and Computer Science. University of Leicester. Leicester, UK.
2	Outline of talk. Description of Higher Order Symmetric Variable Stepsize Methods (HOSVS). Analysis of the Adaptive Verlet method. Implementation of HOSVS Methods. Comparison of the proposed HOSVS Methods. Typical applications for HOSVS methods: Arenstorf Orbits/Solar system. Acknowledgements.
3	HOSVS Methods. HOSVS time-reversible integrators are developed using a combination of time-transformation and composition methods. Composition methods typically used are based on the work of Yoshida [1]/McLachlan [2] and others, where a number of steps with a fixed time step ratio are composed to cancel out errors. Time transformation methods considered are based on the application of a Sundman re-parameterisation applied to the Verlet method. It has been observed that a straightforward approach of composing Adaptive Verlet stages restricts the order of the resulting method to four.
4	Composition Methods. By expanding a method as a composition of exponentials with coefficients we can use the Baker-Campbell-Hausdorf (BCH) Theorem to establish constraints on the coefficients that must be satisfied to obtain a higher order method as described in [1] and [2]. The expansion for the composition of two steps is: Given a second order symmetric method Y_D_twith stepsize Dt, we can obtain a higher order method Y^r_D_t of order r and step size Dt by concatenating Y_D_t with coefficients w₁,w₂,…,w_m to get: The derivation of this is based on the following assumptions: Given a system with flows exp(tX), where X denotes a vector field on some space with coordinates z, time t and system initial conditions z⁰ i.e If we can write X=A + B, we have a map f such that: Then we can increase the order to r, say, by composing several stages to get:
5	The Adaptive Verlet Method. The Verlet method for separable Hamiltonian H=T(p)+V(q) is: The Adaptive Verlet method with re-parameterization function dt/dt=G(q,p) is then defined, in [3], as:
6	Backward Error Analysis of the AV method. A Backward Error analysis is based on the idea that we can find a modified vector field for which the numerical method provides an exact solution, a backward error analysis for the AV method is given in [4]. Given a mechanical Hamiltonian system of the form: where q,p2 R^N, M2 R^NxN is a positive definite matrix, and F=- 5_q V for some smooth, real valued potential function V. If we apply to the system a time transformation of the form: The method can be formally written as: where t_n=nt and, This backward error analysis shows that the numerical solution is a superposition of a smooth function with oscillatory terms dependent on the scalar value gⁿ.
7	Analysis of the HOSVS method. The backward error analysis indicates that the numerical method Y, viewed as a phase-space mapping, is dependent on both the step size t and the scalar value gⁿ, and will be denoted Y_t_,gn. It is possible to derive a general equation for Y_t_,gn as follows: If Y_t is the flow map for the Adaptive Verlet method, at a step size of t, given by, We can write the steps for re-parameterization scalar gⁿ and gⁿ⁺¹ as for some a₁,a₂,…,b₁,b₂,….
8	Analysis of the HOSVS method. We can use our knowledge of the fact that Y is time reversible, from [3], to derive the relationships between a₁,a₂,…, and b₁,b₂,…, Then, using the BCH formula together with the fact that all coefficients of t must be zero for the identity mapping, we have that, (1) requires that a₁-b₁=0, hence b₁=a₁. Thus. (1) also requires that a₂+b₂=0 hence b₂=-a₂, hence, etc. Giving the general equation,
9	Odd-Even Composition. Solution 1: The Odd-Even composition method (OEC). The first solution is to compose the odd and even steps to create a new method on which to base a composition scheme i.e. Using the BCH equation with the expansion for Y_t_,gn we get: Since only odd order terms in t exist in the combined method it is suitable for use in any Yoshida type composition scheme.
10	Cancellation of Even Order Terms. Solution 2: Cancellation of even order terms (CEOT method). By calculating the coefficients recursively, using Yoshida’s approach [1], the new constraints allow for coefficients which cancel the even order terms introduced by the AV method. To identify the constraints, If we generate a method of high order, say r^th order composed of 2m+1 steps of Y, using the BCH formula we find the following, where f₄,f_5b,f₅ are homogeneous polynomials of degrees 4,5,5 respectively. The conditions for all methods with order greater than 4 are:
11	Cancellation of Even Order Terms. The constraints containing even order terms and f₄, f_5b are in addition to those required for a fixed step size scheme. Coefficients can then be calculated using the Newton-Raphson method, one possible set for a 6^th order method is given below:
12	Pre/post processing. Solution 3: Pre/post processing (Mixed Order or MO method). A sequence of steps of the AV method can be rearranged as follows. This produces a Mixed Order method, with q,p evolving with r^th order but the time re-parameterisation evolving with 2^nd order.
13	Comparison of Methods. The MO mixed-order method is high order in q,p but second order in the re-scaling factor g and this can lead to some problems. For high eccentricity Kepler models the error constant is around 80 times worse than the OEC and CEOT methods. For 6^th order methods the CEOT method is more efficient than the OEC method, however for orders greater than 8^th the CEOT method requires more coefficients than the OEC method which then becomes the more efficient. A graph of the two methods follows:
14	Arenstorf Orbit model. The Arenstorf orbit, also known as the restricted three body problem, was used to test the OEC and CEOT methods. It consists of a pair of bodies, m₁ and m₂ , rotating around the z-axis with a third body, of small mass, orbiting the system in the x-y plane. The system units are such that sum of masses, G and angular velocity are one, then, This system has the feature of close approaches between the orbiting body and the smaller of the rotating pair of bodies. In this setting large improvements in integrator efficiency can be obtained by using a variable stepsize or time re-parameterisation based on the distance between the bodies of interest.
15	Arenstorf Orbit simulation.
16	Solar System model. The Solar system is a familiar model and a useful test bed for evaluating integrators. The eccentricity for Mercury is sufficient to produce improvements in efficiency of around 1.8 when using a variable stepsize method based on the Sun-Mercury distance. The function used was, The MO method lends itself to this model as a second order approximation to the re-parameterisation is sufficiently accurate. A graph of the efficiency of the method compared to the fixed step method is given below. Simulations of over 1 Billion years have been done with this method. The Hamiltonian equation used, where m_i is the mass of the i^th planet and G is the gravitational constant is,
17	Solar System Simulation.
18	Conclusions. HOSVS methods provide a useful method of combining the efficiency gains of both higher order and variable stepsize methods. The CEOT method is the most efficient where 6^th-8^th order integrators, which accurately calculate the re-parameterisation function, are required. The OEC method is useful for very high order integrators or where existing coefficients are to be used, which accurately calculate the re-parameterisation function. The MO method is best suited, and the most efficient, where accurate determination of the re-parameterisation function is not important, such as time transformations based on the eccentricity of the orbit of Mercury.
19	Acknowledgements. [1] Yoshida, H.,Construction of higher order symplectic integrators, Physics Letters A, 150, No. 5,6,7, pp 262-268 , 1990. [2] McLachlan, R.,On the numerical integration of ordinary differential equations by symmetric composition methods, SIAM Journal Sci. Comput., 16, No. 1, pp 151-168, 1995. [3] Huang, W. and Leimkuhler, B.,The Adaptive Verlet Method, SIAM Journal Sci. Comput. 18, pp 239-256, 1997. [4] Cirilli, S., Hairer, E., Leimkuhler, B.,Asymptotic Error Analysis of the Adaptive Verlet Method, BIT, 39, pp 25-33, 1999. EPRSC-Student funding.