Leap frog algorith for solving two body problem

Question

I made a program to solve the two body problem with the leap frog algorithm but the output diagram is wrong. I solved the problem by using an other technique (runge kutta, second order or midpoint) with the same initial conditions, and the result was correct. So the problem is not in the initial conditions or in the setting of the matrices, but in how I used the technique. Here is the code:

import numpy  as np
import random as rn 
import math 
from matplotlib import pyplot as plt

#A simple two body problem with Newton's law

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~ ~ ~ ~ ~ ~ FUNCTIONS ~ ~ ~ ~ ~ ~

#This is the acceleration: 

def a(j,i,2): 
#j is the index of a specific particle
#i is the index of the vector
    s = 0
    for k in range(2):
        if (k!=j):
            a = -G*m[k]*((r[j,i]-r[k,i])/pow(np.linalg.norm(r[j,i]-r[k,i]),3))
            s += a
    return s

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~ ~ ~ ~ ~ ~ MAIN BODY ~ ~ ~ ~ ~ ~

Dt = 300
h  = 0.01
N  = int(Dt/h)
G  = 1                  #units: Newton's G is equal to 1

m  = np.ones (2, float) #both masses equal to 1
t  = np.zeros(N, float)
r  = np.zeros((2,N,2), float) #position 
u  = np.zeros((2,N,2), float) #velocity

#first dimension is the two particles and the second
#is the value of x or y, depending on the third dimension
#For example: 
#  r[0,43,0]: First particle, 43th step, x dimension
#  r[1,43,0]: Second particle,43th step, x dimension
#  r[1,43,1]: Second particle,43th step, y dimension
#It is the same for the velocity u

#Initial conditions
r[0,0] = [ 1.0, 1.0]  
r[1,0] = [-1.0,-1.0]
u[0,0] = [-0.5, 0.0]
u[1,0] = [ 0.5, 0.0]

for i in range(1,N,1):
    t[i] = t[i-1] + h
    for j in range(2):
        r[j,i] = r[j,i-1] + h*u[j,i-1]+(0.5*h**2)*a(j,i-1,2)
        u[j,i] = u[j,i-1] + 0.5*h*(a(j,i-1,2)+a(j,i,2)) 

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~ ~ ~ ~ ~ ~ ~ PLOTS ~ ~ ~ ~ ~ ~ ~

plt.plot(r[0,:,0], r[0,:,1], 'black')
plt.plot(r[1,:,0], r[1,:,1],   'red')
plt.show()

Here is the equations of the leap frog technique

score 1 · Answer 1 · answered Dec 09 '20 at 14:16

In the Verlet loop (velocity, not leapfrog)

for i in range(1,N,1):
    t[i] = t[i-1] + h
    for j in range(2):
        r[j,i] = r[j,i-1] + h*u[j,i-1]+(0.5*h**2)*a(j,i-1,2)
        u[j,i] = u[j,i-1] + 0.5*h*(a(j,i-1,2)+a(j,i,2))

you are using the still-to-compute position r[1,i] in the computation of u[0,i]. This will of course give a wrong result. You need to make two loops, all positions first, then velocities

for i in range(1,N,1):
    t[i] = t[i-1] + h
    for j in range(2):
        r[j,i] = r[j,i-1] + h*u[j,i-1]+(0.5*h**2)*a(j,i-1,2)
    for j in range(2):
        u[j,i] = u[j,i-1] + 0.5*h*(a(j,i-1,2)+a(j,i,2))

Leap frog algorith for solving two body problem

1 Answers1