implemented symmetric high-order symplectic decomposition by Yoshida

examples/test_symp_mo_1p_1d.py

+"""  High-order symplectic integrators """
+import numpy as np
+import matplotlib.pyplot as plt
+from rotagaporp.symplectic import HighOrderSymplectic, Symplectic
+from rotagaporp.systems import Morse
+
+syspar = {'q0': 3.0, 'p0': 0.0}
+system = Morse(**syspar)
+
+params = {'tstart': 0.0, 'tend': 10.0, 'tstep': 0.1, 'tqdm': True}
+prop = HighOrderSymplectic(order=4, syst=system, **params)
+# prop = Symplectic(syst=system, **params)
+prop.propagate()
+prop.analyse()
+(ti, qt, pt) = prop.get_trajectory()
+print('energy drift of symplectic propagator', np.max(np.fabs(prop.er)))
+
+q_an, p_an = system.solution(times=ti)
+en0 = system.energy(system.q0, system.p0)
+evecf = np.vectorize(system.energy)
+err_an = np.array((evecf(q_an, p_an)-en0)/en0, dtype=float)
+print('energy drift of analytic solution', np.max(np.fabs(err_an)))
+
+fig = plt.figure()
+plot = fig.add_subplot(311)
+# plt.plot(ti, qt, label='q symplectic')
+# plt.plot(ti, q_an, label='q analytic')
+plt.plot(ti, abs(qt-q_an), label='q diff')
+plot.set_ylabel(r'$q(t)$')
+plt.legend()
+
+plot = fig.add_subplot(312)
+# plt.plot(ti, pt, label='p symplectic')
+# plt.plot(ti, p_an, label='p analytic')
+plt.plot(ti, abs(pt-p_an), label='p diff')
+plot.set_ylabel(r'$p(t)$')
+plt.legend()
+
+plot = fig.add_subplot(313)
+plt.plot(ti, np.maximum.accumulate(abs(prop.er)), label='energy drift')
+plot.set_ylabel(r'$\Delta E$')
+plot.set_xlabel('time')
+plt.legend()
+
+plt.show()

rotagaporp/symplectic.py

 """ Symplectic propagator for classical particle dynamics """
+import sympy as sp
 from propagators import Propagator
 
 class Symplectic(Propagator):
@@ -31,3 +32,78 @@ class Symplectic(Propagator):
         self.force_saved = self.syst.force(qt)
         pt = p_half + 0.5*self.tstep*self.force_saved
         return (qt, pt)
+
+
+class HighOrderSymplectic(Propagator):
+    """ High-order symplectic propagator for classical particle dynamics """
+    name = 'high order symplectic'
+
+    def __init__(self, order=None, **kwargs):
+        super().__init__(**kwargs)
+        assert order > 0, 'order must be positive'
+        assert order%2 == 0, 'order must be even integer'
+        self.order = order
+        self.set_expression_1()
+
+    def set_expression_1(self):
+        """ First version, H. Yoshida, Phys. Lett. A 150, 262 (1990) """
+        A = sp.diff(self.syst.potential_energy, self.syst.q)
+        B = sp.diff(self.syst.kinetic_energy, self.syst.p)
+        self.qt = self.syst.q
+        self.pt = self.syst.p
+        # the c and d coefficients in the cited paper are interchanged
+        for ci, di in reversed(list(zip(*self.get_coefficients()))):
+            self.qt += self.tstep*di*B.subs(self.syst.p, self.pt)
+            self.pt -= self.tstep*ci*A.subs(self.syst.q, self.qt)
+
+    def get_coefficients(self):
+        """ Decomposition coefficients prod_i exp(ci*tau*A)*exp(di*tau*B)
+            H. Yoshida, Phys. Lett. A 150, 262 (1990) """
+        if self.order > 2:
+            A = sp.Symbol('A', commutative=False)
+            B = sp.Symbol('B', commutative=False)
+            tau = sp.Symbol('tau', real=True)
+            S = self.get_decomposition(A, B, tau)
+            ccoef = [i.args[0].args[0] for i in S.args if A in i.args[0].args]
+            dcoef = [i.args[0].args[0] for i in S.args if B in i.args[0].args]
+            dcoef.append(0)
+        else:
+            ccoef = [0.5, 0.5]
+            dcoef = [1, 0]
+        return ccoef, dcoef
+
+    def set_expression_2(self):
+        """ Second version, H. Yoshida, Phys. Lett. A 150, 262 (1990) """
+        q = self.syst.q
+        p = self.syst.p
+        DV = -sp.diff(self.syst.potential_energy, q)
+        DT = sp.diff(self.syst.kinetic_energy, p)
+        A = sp.Symbol('A', commutative=False)
+        B = sp.Symbol('B', commutative=False)
+        tau = sp.Symbol('tau', real=True)
+        self.qt = q
+        self.pt = p
+        for s in self.get_decomposition(A, B, tau).args:
+            if A in s.args[0].args:
+                self.qt = self.qt.subs({q: q+s.args[0].subs({A: DT})})
+                self.pt = self.pt.subs({q: q+s.args[0].subs({A: DT})})
+            else:
+                self.qt = self.qt.subs({p: p+s.args[0].subs({B: DV})})
+                self.pt = self.pt.subs({p: p+s.args[0].subs({B: DV})})
+        self.qt = self.qt.xreplace({tau: self.tstep})
+        self.pt = self.pt.xreplace({tau: self.tstep})
+
+    def get_decomposition(self, A, B, tau):
+        """ Symmetric decomposition of the propagator exp(tau*(A+B))
+            H. Yoshida, Phys. Lett. A 150, 262 (1990) """
+        S = sp.exp(tau*A/2)*sp.exp(tau*B)*sp.exp(tau*A/2)
+        for n in range(1, self.order//2):
+            z0 = -2**(1/(2*n+1))/(2-2**(1/(2*n+1)))
+            z1 = 1/(2-2**(1/(2*n+1)))
+            S = S.subs(tau, z1*tau)*S.subs(tau, z0*tau)*S.subs(tau, z1*tau)
+        return S
+
+    def step_nopbc(self, time, q0, p0):
+        """ Time stepping """
+        repl = {self.syst.q: q0, self.syst.p: p0}
+        return self.qt.xreplace(repl), self.pt.xreplace(repl)