time-symmetric (time reversible) propagator and example

examples/test_cheb_mo_1p_1d_symmetric.py

+""" One particle in Morse potential - Chebyshev and symmetrized Chebyshev """
+import numpy as np
+import matplotlib.pyplot as plt
+from rotagaporp.chebyshev import Chebyshev, SymmetricChebyshev
+from rotagaporp.symplectic import Symplectic
+from rotagaporp.systems import Morse
+from generic import spfloat
+
+prec = 30
+
+params = {'tstart': spfloat(0., prec), 'tend': spfloat(100., prec),
+          'tstep': spfloat(0.02, prec), 'nterms': 4, 'tqdm': True}
+syspar = {'q0': spfloat(3., prec), 'p0': spfloat(0., prec),
+           'fpprec': 'multi', 'prec': prec}
+system = Morse(**syspar)
+
+# velocity Verlet
+params_verlet = params.copy()
+del params_verlet['nterms']
+vv_prop = Symplectic(syst=system, **params_verlet)
+vv_prop.propagate()
+vv_prop.analyse()
+
+params['fpprec'] = 'multi'
+params['prec'] = prec
+
+# one-step propagator
+ch_prop = Chebyshev(syst=system, **params)
+ch_prop.propagate()
+time_ch, q_ch, p_ch = ch_prop.get_trajectory()
+ch_prop.analyse()
+print('energy drift of one-step prop', np.max(np.fabs(ch_prop.er)))
+
+# symmetrized propagator
+params['signed'] = False
+ts_prop = SymmetricChebyshev(syst=system, **params)
+ts_prop.propagate()
+time_ts, q_ts, p_ts = ts_prop.get_trajectory()
+ts_prop.analyse()
+print('energy drift of symmetric prop', np.max(np.fabs(ts_prop.er)))
+
+assert len(time_ts) == len(time_ch)
+assert np.allclose(time_ts, time_ch)
+q_an, p_an = system.solution(times=time_ts)
+
+fig = plt.figure()
+plot1 = fig.add_subplot(311)
+plot1.plot(time_ch, np.maximum.accumulate(abs(ch_prop.er)), label='Chebyshev')
+plot1.plot(time_ts, np.maximum.accumulate(abs(ts_prop.er)), label='Symmetric Chebyshev')
+plot1.plot(time_ts, np.maximum.accumulate(abs(vv_prop.er)), label='Velocity Verlet')
+plot1.set_ylabel('maximum energy drift')
+plot1.set_yscale('log')
+plot1.set_xscale('log')
+plot1.legend()
+
+plot2 = fig.add_subplot(312)
+plot2.plot(time_ts, q_ch, label='Chebyshev')
+plot2.plot(time_ts, q_ts, label='Symmetric Chebyshev')
+plot2.plot(time_ts, q_an, label='Analytic solution')
+plot2.set_xlabel(r'$t$')
+plot2.set_ylabel(r'$q(t)$')
+plot2.set_xscale('log')
+plot2.legend()
+
+plot3 = fig.add_subplot(313)
+plot3.plot(time_ts, p_ch, label='Chebyshev')
+plot3.plot(time_ts, p_ts, label='Symmetric Chebyshev')
+plot3.plot(time_ts, p_an, label='Analytic solution')
+plot3.set_xlabel(r'$t$')
+plot3.set_ylabel(r'$p(t)$')
+plot3.set_xscale('log')
+plot3.legend()
+
+plt.show()

rotagaporp/chebyshev.py

@@ -5,6 +5,7 @@ import sympy as sp
 from scipy.special import jn, iv
 from propagators import PolynomialPropagator
 from propagators import TwoStepPolynomialPropagator
+from propagators import SymmetricPolynomialPropagator
 
 
 def cheb_coef(bf, nterms, alpha):
@@ -97,3 +98,24 @@ class TwoStepChebyshev(TwoStepPolynomialPropagator):
         self.coeffs = 2.*cheb_coef(bf, self.nterms, alpha)[self.parity::2]
         self.set_polynomial(system=self.syst, scale=2./self.delta, sign=sign,
                             recursion='chebyshev', parity=self.parity)
+
+
+class SymmetricChebyshev(SymmetricPolynomialPropagator):
+    """ Symmetric Chebyshev propagator for classical particle dynamics """
+    name = 'symmetric chebyshev'
+    def __init__(self, delta=1., nterms=None, signed=False, **kwargs):
+        super().__init__(**kwargs)
+
+        assert delta > 0
+        self.delta = delta
+        alpha = self.delta*self.tstep/2.
+        self.nterms = int(np.ceil(alpha)) if nterms is None else nterms
+        self.efficiency = alpha/self.nterms
+        bf = besself(self.fpprec, signed, prec=self.prec)
+        assert isinstance(signed, bool)
+        sign = -1 if signed else 1
+        assert self.substep is None, 'substep not implemented'
+        self.coeffs = cheb_coef(bf, self.nterms, alpha)
+        self.reverse_coeffs = cheb_coef(bf, self.nterms, -alpha)
+        self.set_polynomial(system=self.syst, scale=2./self.delta, sign=sign,
+                            recursion='chebyshev')

rotagaporp/propagators.py

@@ -126,12 +126,10 @@ class PolynomialPropagator(Propagator):
         pt = np.dot(self.coeffs, pvec)
         return (qt, pt)
 
-    def step_pbc(self, time, qval, pval):
+    def step_pbc(self, *args, **kwargs):
         """ perform a time step with periodic boundary conditions """
-        qvec, pvec = self.iln.get_val_float(qval, pval)
-        qt = np.dot(self.coeffs, qvec)
-        pt = np.dot(self.coeffs, pvec)
-        return (self.syst.wrap_q(qt), pt)
+        wrap = lambda x, y: (self.syst.wrap_q(x), y)
+        return wrap(*self.step_nopbc(*args, **kwargs))
 
     def load(self, filename='propag.pkl'):
         """ deserialize the propagator object and extract the iln object """
@@ -181,11 +179,27 @@ class TwoStepPolynomialPropagator(PolynomialPropagator):
         self.p_1 = pval
         return (qt, pt)
 
-    def step_pbc(self, time, qval, pval):
-        """ perform a time step with periodic boundary conditions """
+
+class SymmetricPolynomialPropagator(PolynomialPropagator):
+    """ Self-starting one-step time-symmetric polynomial propagators """
+
+    name = 'symmetric polynomial propagator'
+    iln = None
+    coeffs = None
+    reverse_coeffs = None
+
+    def __init__(self, **kwargs):
+        """ """
+        super().__init__(**kwargs)
+
+    def step_nopbc(self, time, qval, pval):
+        """ perform a time step with initial coordinate and momentum """
         qvec, pvec = self.iln.get_val_float(qval, pval)
-        qt = np.dot(self.coeffs, qvec) + self.prefact*self.q_1
-        pt = np.dot(self.coeffs, pvec) + self.prefact*self.p_1
-        self.q_1 = qval
-        self.p_1 = pval
-        return (self.syst.wrap_q(qt), pt)
+        q_for = np.dot(self.coeffs, qvec)
+        p_for = np.dot(self.coeffs, pvec)
+        qvec, pvec = self.iln.get_val_float(q_for, p_for)
+        q_rev = np.dot(self.reverse_coeffs, qvec)
+        p_rev = np.dot(self.reverse_coeffs, pvec)
+        q_res = qval + q_for - q_rev
+        p_res = pval + p_for - p_rev
+        return (q_res, p_res)