Automatic code generator for higher order integrators

Abstract

Some explicit algorithms for higher order symplectic integration of a large class of Hamilton’s equations have recently been discussed by Mushtaq et al. Here we present a Python program for automatic numerical implementation of these algorithms for a given Hamiltonian, both for double precision and multiprecision computations. We provide examples of how to use this program, and illustrate behavior of both the code generator and the generated solver module(s). Program summaryProgram title: HOMsPy: Higher Order (Symplectic) Methods in PythonCatalogue identifier: AESD_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AESD_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 19423No. of bytes in distributed program, including test data, etc.: 1970283Distribution format: tar.gzProgramming language: Python 2.7.Computer: PCs or higher performance computers.Operating system: Linux, MacOS, MSWindows.RAM: Kilobytes to a several gigabytes (problem dependent).Classification: 4.3, 5.External routines: SymPy library [1] for generating the code. NumPy library [2], and optionally mpmath [3] library for running the generated code. The matplotlib [4] library for plotting results.Nature of problem:We have developed algorithms [5] for numerical solution of Hamilton’s equations. (1)q̇a=∂H(q,p)/∂paa,ṗa=−∂H(q,p)/∂qa,a=1…N for Hamiltonians of the form (2)H(q,p)=T(p)+V(q)=(1/2)pTMp+V(q), with M a symmetric positive definite matrix. The algorithms preserve the symplectic property of the time evolution exactly, and are of orders τN (for 2≤N≤8) in the timestep τ. Although explicit, the algorithms are time-consuming and error-prone to implement numerically by hand, in particular for larger N.Solution method:We use computer algebra to perform all analytic calculations required for a specific model, and to generate the Python code for numerical solution of this model, including example programs using that code.Restrictions:In our implementation the mass matrix is assumed to be equal to the unit matrix, and V(q) must be sufficiently differentiable.Running time:Subseconds to eons (problem dependent). See discussion in the main article.