implement the leapfrog and stormer_verlet in jobs/src/integrators.py

jobs/src/integrators.py

@@ -20,3 +20,53 @@ def verlet(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, dt: float
     p = p + 1 / 2 * F(q) * dt
 
     return p, q
+
+
+def leapfrog(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, dt: float):
+    '''
+        # using leapfrog integration. https://en.wikipedia.org/wiki/Leapfrog_integration
+        The basic idea of the Leapfrog method is to update the position
+        and momentum of a particle or system at half-integer time steps.
+        The algorithm proceeds in discrete steps, and the update scheme looks like this:
+        1. Update Velocity at Half Step:
+        2. Update Position:
+        3. Update Velocity at Full Step:
+    '''
+    p = p0
+    q = q0
+
+    # 1. Update Momentum at Half Step:
+    p += dt / 2 * F(q)
+
+    # 2. Update Position:
+    q += p / m * dt
+
+    # 3. Update Momentum at Full Step:
+    p += dt / 2 * F(q)
+
+    # print("-------------------")
+    # print(f'accelerations={accelerations}')
+    # print(f'velocity={velocity}')
+    return p, q
+
+
+
+def stormer_verlet(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, dt: float):
+    """
+      stormer_verlet integrator for one time step
+      # https://www.physics.udel.edu/~bnikolic/teaching/phys660/numerical_ode/node5.html
+      # https://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet
+      # https://www.algorithm-archive.org/contents/verlet_integration/verlet_integration.html
+    """
+    p = p0
+    q = q0
+
+    # 1. Update Position at Half Step:
+    acc = F(q) / m
+    q += p / m * dt + acc * (dt ** 2 * 0.5)
+
+    # 2. Update Momentum:
+    p += acc * (dt * 0.5) * m + F(q) * (dt * 0.5)
+
+    return p, q
+