implement leapfrog and stormer_verlet
5 unresolved threads
5 unresolved threads
implement leapfrog and stormer_verlet
Merge request reports
Activity
requested review from @nguyed99
assigned to @jim92
35 ''' 36 p = p0 37 q = q0 38 39 # 1. Update Momentum at Half Step: 40 p = p + dt / 2 * F(q) 41 42 # 2. Update Position: 43 q = q + p / m * dt 44 45 # 3. Update Momentum at Full Step: 46 p = p + dt / 2 * F(q) 47 48 # print("-------------------") 49 # print(f'accelerations={accelerations}') 50 # print(f'velocity={velocity}') 20 20 p = p + 1 / 2 * F(q) * dt 21 21 22 22 return p, q 23 24 25 26 def leapfrog(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, dt: float): 27 ''' 28 # using leapfrog integration. https://en.wikipedia.org/wiki/Leapfrog_integration 29 The basic idea of the Leapfrog method is to update the position 30 and momentum of a particle or system at half-integer time steps. 45 # 3. Update Momentum at Full Step: 46 p = p + dt / 2 * F(q) 47 48 # print("-------------------") 49 # print(f'accelerations={accelerations}') 50 # print(f'velocity={velocity}') 51 return p, q 52 53 54 55 def stormer_verlet(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, dt: float): 56 """ 57 stormer_verlet integrator for one time step 58 # https://www.physics.udel.edu/~bnikolic/teaching/phys660/numerical_ode/node5.html 59 # https://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet 60 # https://www.algorithm-archive.org/contents/verlet_integration/verlet_integration.html 55 def stormer_verlet(F: Callable, q0: np.ndarray, p0: np.ndarray, m: np.ndarray, dt: float): 56 """ 57 stormer_verlet integrator for one time step 58 # https://www.physics.udel.edu/~bnikolic/teaching/phys660/numerical_ode/node5.html 59 # https://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet 60 # https://www.algorithm-archive.org/contents/verlet_integration/verlet_integration.html 61 """ 62 p = p0 63 q = q0 64 65 # 1. Update Position at Half Step: 66 acc = F(q) / m 67 q = q + p / m * dt + acc * (dt ** 2 * 0.5) 68 69 # 2. Update Momentum: 70 p = p + acc * (dt * 0.5) * m + F(q) * (dt * 0.5) 31 The algorithm proceeds in discrete steps, and the update scheme looks like this: 32 1. Update Velocity at Half Step: 33 2. Update Position: 34 3. Update Velocity at Full Step: 35 ''' 36 p = p0 37 q = q0 38 39 # 1. Update Momentum at Half Step: 40 p = p + dt / 2 * F(q) 41 42 # 2. Update Position: 43 q = q + p / m * dt 44 45 # 3. Update Momentum at Full Step: 46 p = p + dt / 2 * F(q)
Please register or sign in to reply