Skip to content
Snippets Groups Projects
Commit 41e489af authored by JayM0826's avatar JayM0826
Browse files

fix bug for reverse unit test.

All tests in local passed
parent 907a611f
No related branches found
No related tags found
1 merge request!7leapfrog and stormer_verlet
Pipeline #59214 passed
...@@ -22,6 +22,7 @@ def verlet(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, dt: float ...@@ -22,6 +22,7 @@ def verlet(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, dt: float
return p, q return p, q
def leapfrog(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, dt: float): 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 # using leapfrog integration. https://en.wikipedia.org/wiki/Leapfrog_integration
...@@ -36,13 +37,13 @@ def leapfrog(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, dt: flo ...@@ -36,13 +37,13 @@ def leapfrog(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, dt: flo
q = q0 q = q0
# 1. Update Momentum at Half Step: # 1. Update Momentum at Half Step:
p += dt / 2 * F(q) p = p + dt / 2 * F(q)
# 2. Update Position: # 2. Update Position:
q += p / m * dt q = q + p / m * dt
# 3. Update Momentum at Full Step: # 3. Update Momentum at Full Step:
p += dt / 2 * F(q) p = p + dt / 2 * F(q)
# print("-------------------") # print("-------------------")
# print(f'accelerations={accelerations}') # print(f'accelerations={accelerations}')
...@@ -63,10 +64,10 @@ def stormer_verlet(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, d ...@@ -63,10 +64,10 @@ def stormer_verlet(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, d
# 1. Update Position at Half Step: # 1. Update Position at Half Step:
acc = F(q) / m acc = F(q) / m
q += p / m * dt + acc * (dt ** 2 * 0.5) q = q + p / m * dt + acc * (dt ** 2 * 0.5)
# 2. Update Momentum: # 2. Update Momentum:
p += acc * (dt * 0.5) * m + F(q) * (dt * 0.5) p = p + acc * (dt * 0.5) * m + F(q) * (dt * 0.5)
return p, q return p, q
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment