Symplectic Integration

Synopsis

#include <boost/math/quadrature/symmetric.hpp>
namespace boost{ namespace math{ namespace quadrature {

enum class available_methods
{
    Y6,
    Y4,
    Y2,
    SRKNB6,
    SRKNB11
};

template <typename RandomAccessContainer, typename RealType, class Func, class Policy>
std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian(RandomAccessContainer p0,
                                                                                                         RandomAccessContainer q0,
                                                                                                         const RealType dt,
                                                                                                         const unsigned steps,
                                                                                                         Func dHdp,
                                                                                                         Func dHdq,
                                                                                                         available_methods method,
                                                                                                         const Policy& pol)

template <typename RandomAccessContainer, typename RealType, class Func>
std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian(RandomAccessContainer p0,
                                                                                                         RandomAccessContainer q0,
                                                                                                         const RealType dt,
                                                                                                         const unsigned steps,
                                                                                                         Func dHdp,
                                                                                                         Func dHdq,
                                                                                                         available_methods method)

template <typename RandomAccessContainer, typename RealType, class Func>
std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian(RandomAccessContainer p0,
                                                                                                         RandomAccessContainer q0,
                                                                                                         const RealType dt,
                                                                                                         const unsigned steps,
                                                                                                         Func dHdp,
                                                                                                         Func dHdq)

template <typename RandomAccessContainer, typename RealType, class Func, class Policy>
BOOST_MATH_EXPORT std::tuple<std::vector<RealType>, std::vector<RealType>, std::vector<RealType>, std::vector<RealType> > integrate_hamiltonian_adaptive(RandomAccessContainer p0,
                                                                                                                                                         RandomAccessContainer q0,
                                                                                                                                                         const std::pair<RealType, RealType> timeInterval,
                                                                                                                                                         const Func dHdp,
                                                                                                                                                         const Func dHdq,
                                                                                                                                                         const RealType atol,
                                                                                                                                                         const RealType rtol,
                                                                                                                                                         const available_methods& method,
                                                                                                                                                         const Policy pol)
template <typename RandomAccessContainer, typename RealType, class Func>
BOOST_MATH_EXPORT std::tuple<std::vector<RealType>, std::vector<RealType>, std::vector<RealType>, std::vector<RealType> > integrate_hamiltonian_adaptive(RandomAccessContainer p0,
                                                                                                                                                        RandomAccessContainer q0,
                                                                                                                                                        const std::pair<RealType, RealType> timeInterval,
                                                                                                                                                        const Func dHdp,
                                                                                                                                                        const Func dHdq,
                                                                                                                                                        const RealType atol,
                                                                                                                                                        const RealType rtol,
                                                                                                                                                        const available_methods& method)

template <typename RandomAccessContainer, typename RealType, class Func>
BOOST_MATH_EXPORT std::tuple<std::vector<RealType>, std::vector<RealType>, std::vector<RealType>, std::vector<RealType> > integrate_hamiltonian_adaptive(RandomAccessContainer p0,
                                                                                                                                                        RandomAccessContainer q0,
                                                                                                                                                        const std::pair<RealType, RealType> timeInterval,
                                                                                                                                                        const Func dHdp,
                                                                                                                                                        const Func dHdq,
                                                                                                                                                        const RealType atol,
                                                                                                                                                        const RealType rtol)

template <typename RealType, class Func, class Policy>
BOOST_MATH_EXPORT boost::math::interpolators::cubic_hermite<std::vector<RealType> > integrate_hamiltonian_cubic(RealType p0,
                                                                                                                RealType q0,
                                                                                                                const RealType dt,
                                                                                                                const unsigned steps,
                                                                                                                Func dHdp,
                                                                                                                Func dHdq,
                                                                                                                available_methods method,
                                                                                                                const Policy& pol)

template <typename RealType, class Func>
BOOST_MATH_EXPORT boost::math::interpolators::cubic_hermite<std::vector<RealType> > integrate_hamiltonian_cubic(RealType p0,
                                                                                                                RealType q0,
                                                                                                                const RealType dt,
                                                                                                                const unsigned steps,
                                                                                                                Func dHdp,
                                                                                                                Func dHdq,
                                                                                                                available_methods method)

template <typename RealType, class Func>
BOOST_MATH_EXPORT boost::math::interpolators::cubic_hermite<std::vector<RealType> > integrate_hamiltonian_cubic(RealType p0,
                                                                                                                RealType q0,
                                                                                                                const RealType dt,
                                                                                                                const unsigned steps,
                                                                                                                Func dHdp,
                                                                                                                Func dHdq)
}}} // namespaces

The above functions calculate the phase space trajectory for a given Hamiltonian. The functional integrate_hamiltonian integrates a multi-dimensional system given a timestep and number of steps and returns the phase-space trajectory of (p, q) at each time. Then a Cubic Hermite spline can be used to interpolate to any p using the derivative dp/dt = -dH/dq. This leads to the next functional, integrate_hamiltonian_cubic, which integrates a 1 dimensional system given a timestep and number of steps and returns a Cubic Hermite spline which interpolates p to any time. Lastly, the functional integrate_hamiltonian_adaptive integrates a multi-dimensional system to a given absolute and relative tolerance on a given time interval, (tmin, tmax), using an integrating controller to adjust the timestep. This function returns the tuple (t, p, q, error).

The trajectories are calculated using symplectic integration which preserves the energy of a system. Even higher order traditional numerical integration algorithms will gain or lose energy at long times. The functional implements the methods in Yoshida's landmark paper and methods by Blanes and Moan. We assume that the Hamiltonian is separable so that it can be written as H = T(p) + V(q).

We now give an example for a simple harmonic oscillator with the Hamiltonian

For simplicity we will assume k = m = 1. Then the partial derivatives of the Hamiltonian with respect to p and q are

double dHdp(double p)
{
    return p;
}

double dHdq(double q)
{
    return q;
}

Note that the functional can be readily generalized to multiple coordinates by changing the signature of dHdp to

std::vector<double> dHdp(std::vector<double> p)
{
    // calculate the partial derivatives with respect to each p_i
}

We then define the timestep size and number of steps of the algorithm to go from t=0 to t=100

RealType dt = 0.05;
RealType t_end = 100;
unsigned int steps = t_end / dt;

Lastly, we define the initial conditions so that the oscillator starts from rest at x=1

RealType q0 = 1;
RealType p0 = 0;

We now evolve the system using the following

std::vector<RealType> p;
std::vector<RealType> q;

std::tie(p, q) = boost::math::quadrature::integrate_hamiltonian(p0, q0, dt, steps, dHdp, dHdq, available_methods::Y6);

The ouput vectors q, p are the position and momentum of the system at each time. In higher dimensions the output will be a vector of vectors.

The last argument that we pass to integrate_hamiltonian, available_methods::Y6 here, sets the integration method to use. The value Y6 stands for Yoshida's 6th order integrator. Yoshida's 2nd and 4th order methods are also available by passing the value Y2, or Y4 respectively. The 4th order and 6th order symplectic Runge-Kutta-Nystrom methods from Table 3 of Blanes and Moan are available using the values SRKNB6 and SRKNB11.

The final Policy argument is optional and can be used to control the behaviour of the function: how it handles errors, what level of precision to use etc. Refer to the policy documentation for more details.

References:

Yoshida, Haruo. [`Construction of higher order symplectic integrators], Physics Letters A, 150.5-7 (1990): 262-268.

Blanes, Sergio, and Per Christian Moan. [Practical symplectic partitioned runge–kutta and runge–kutta–nyström methods], Journal of Computational and Applied Mathematics 142.2 (2002): 313-330.

Symplectic Integration

Synopsis

Description