I have been attempting to get a two body problem to work, which in the future should be used for more planets, but it is not working and the plot I am supposed to get is circular but I am receiving a straight line for the two body systems. Does anyone know how I can fix this and get the correct plot?

This is the code I use:

import numpy as np
import matplotlib.pyplot as plt

aEarth = 1 #semi-major axis (AU)
eEarth = 0.0167 #eccentricity (no units)

aMercury = 0.387098 #semi-major axis (AU)
eMercury = 0.205635 #eccentricity (no units)
    
Msun = 1 #Mass of Sun (Solar Mass)
Mearth = 3.0024584*10**(-6) #Mass of Earth (Solar Mass)
Mmercury = 1.65956463*10**(-7) #Mass of Mercury (Solar Mass)
Mes = (Msun   Mearth)
Mms = (Msun   Mmercury)
G = 1 #Gravitational Constant 

apocentreEarth = (aEarth*(1   eEarth))
apocentreMercury = (aMercury*(1   eMercury))

vapocentreEarth = (np.sqrt((G*(Mearth Msun)/aEarth)*((1-eEarth)/(1 eEarth))))
vapocentreMercury = (np.sqrt((G*(Mmercury Msun)/aMercury)*(1-eMercury/1 eMercury)))

xEarth = apocentreEarth
yEarth = 0.0
zEarth = 0.0

vxEarth = 0.0
vyEarth =(vapocentreEarth)
vzEarth = 0.0

xMercury = apocentreMercury
yMercury = 0.0
zMercury = 0.0

vxMercury = 0.0
vyMercury =(vapocentreMercury)
vzMercury = 0.0

class CelBody(object):
    # Constants of nature
    def __init__(self, id, name, x0, v0, mass, color, lw):
        # Name of the body (string)
        self.id = id
        self.name = name
        # Mass of the body (solar mass)
        self.M = mass
        # Initial position of the body (au)
        self.x0 = np.asarray(x0, dtype=float)
        # Position (au). Set to initial value.
        self.x = self.x0.copy()
        # Initial velocity of the body (au/s)
        self.v0 = np.asarray(v0, dtype=float)
        # Velocity (au/s). Set to initial value.
        self.v = self.v0.copy()
        self.a = np.zeros([3], dtype=float)
        self.color = color
        self.lw = lw

# All Celestial Objects

t = 0
dt = 0.01

Bodies = [
    CelBody(0, 'Sun', [0, 0, 0], [0, 0, 0], Msun, 'yellow', 10),
    CelBody(1, 'Earth', [xEarth, yEarth, zEarth], [vxEarth, vyEarth, vzEarth], Mearth, 'blue', 3),
    CelBody(2, 'Mercury', [xMercury, yMercury, zMercury], [ vxMercury, vyMercury, vzMercury], Mmercury, 'pink', 3),
    ]

paths = [ [ b.x[:2].copy() ] for b in Bodies]

# loop over ten astronomical years
v = 0
for t in range(0,1000):
    # compute forces/accelerations
    for body in Bodies:
        body.a *= 0
        for other in Bodies:
            # no force on itself
            if (body == other): continue # jump to next loop
            rx = body.x - other.x
            r3 = (np.sqrt(rx[0]**2 rx[1]**2 rx[2]**2))**3
            body.a = -G*other.M*rx/r3

    for n, planet in enumerate(Bodies):
        # use the Forward Euler algorithm 
        planet.a = -G*other.M*rx/r3
        planet.v = planet.v   planet.a*dt
        planet.x = planet.x   planet.v*dt
        paths[n].append( planet.x[:2].copy() )
        #print("s x:Ss v:Ss"%(planet.name,planet.x, planet.v))


plt.figure(figsize=(8,8))
for n, planet in enumerate(Bodies): 
    px, py=np.array(paths[n]).T; 
    plt.plot(px, py, color=planet.color, lw=planet.lw)
plt.show()

CodePudding user response：

There is a problem in these lines:

    planet.v  = planet.v   planet.a*dt
    planet.x  = planet.x   planet.v*dt

They are equivalent to

    planet.v = 2*planet.v   planet.a*dt
    planet.x = 2*planet.x   planet.v*dt

which is not what you want.

Either don't use =, or change those lines to:

    planet.v  = planet.a*dt
    planet.x  = planet.v*dt

There is still a problem: the first line changes planet.v, and then the second used the new planet.v to update planet.x, which is not how the explicit Euler integration should proceed. (There is a comment in the code abou using symplectic Euler, but a comment to this answer says you meant to use forward (i.e. explicit) Euler.)

For this system, a simple fix is to switch the statements:

    planet.x  = planet.v*dt
    planet.v  = planet.a*dt

There may be other issues. If you need more help, please simplify the code as much as possible to create a minimal reproducible example. As it is now, it looks like there are quite a few irrelevant variables declared, two different places where you are computing Euler's method, two places where you assign dt, three bodies defined when you say you are solving the two body problem, etc.

Edit After your update, there are a couple bugs left:

• In the loop that begins for body in Bodies:, in which you compute body.a for each body, the results should be accumulated, so the line that updates body.a should be

              body.a  = -G*other.M*rx/r3
  

  (note the change to =).

• In second inner loop (for n, planet in enumerate(Bodies):), in which you apply the Euler method step, delete the line

          planet.a = -G*other.M*rx/r3
  

  You already computed the accelerations in the previous loop and stored them in planet.a.

If I make just those two changes to the code that is now in the question, the plot shows the circular orbits that you expected.