extended analysis with relative energy difference

classical/propagators.py

@@ -68,9 +68,20 @@ class Propagator:
         pt.shape = (len(pt), -1, 3)
         return ti, qt, pt
 
-    def analyse(self, step=1):
-        """ analyse the trajectory """
+    def analyse(self, step=1, p='err'):
+        """ calculate relative energy error or relative energy difference
+            (p = 1, 2, 3, ..., inf), see refs. https://doi.org/10.2307/2683905
+            https://en.wikipedia.org/wiki/Relative_change_and_difference """
         from itertools import islice
-        qptup = islice(zip(self.sq, self.sp), 0, None, step)
-        self.en = np.array([self.syst.energy(*t) for t in qptup])
-        self.er = np.array((self.en-self.en0)/self.en0, dtype=float)
+        if not hasattr(self, 'en'):
+            qptup = islice(zip(self.sq, self.sp), 0, None, step)
+            self.en = np.array([self.syst.energy(*t) for t in qptup])
+        if p == 'inf':
+            func = lambda x, y: (x-y)/np.maximum(abs(x), abs(y))
+        elif p == 'err':
+            func = lambda x, y: (x-y)/y
+        elif isinstance(p, int) and p > 0:
+            func = lambda x, y: (x-y)/(abs(x)**p+abs(y)**p)**(1/p)
+        else:
+            raise ValueError('invalid value for p: '+repr(p))
+        self.er = np.array(func(self.en, self.en0), dtype=float)