From 907a611ff6a7a99fe5e800030347355352512cca Mon Sep 17 00:00:00 2001
From: JayM0826 <maji91fan@gmail.com>
Date: Mon, 12 Feb 2024 20:21:04 +0100
Subject: [PATCH] implement the leapfrog and stormer_verlet in
 jobs/src/integrators.py

---
 jobs/src/integrators.py | 50 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 50 insertions(+)

diff --git a/jobs/src/integrators.py b/jobs/src/integrators.py
index 35d513a..49d4144 100644
--- a/jobs/src/integrators.py
+++ b/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
+
-- 
GitLab