python class numerical-methods differential-equations orbital-mechanics

Python Euler Method implementation in Two Body Problem not working

I am currently trying to get a two body problem to work, that I can then upgrade to more planets, but it is not working. It is outputting me impossible positions. Does anyone know what is causing that?

This is the code I use:

day = 60*60*24
# Constants
G = 6.67408e-11
dt = 0.1*day
au = 1.496e11
t = 0


class CelBody:

    def __init__(self, id, name, x0, y0, z0, vx0, vy0, vz0, mass, vector, ax0, ay0, az0, totalforcex, totalforcey, totalforcez):
        self.ax0 = ax0
        self.ay0 = ay0
        self.az0 = az0

        self.ax = self.ax0
        self.ay = self.ay0
        self.az = self.az0

        # Constants of nature
        # Universal constant of gravitation
        self.G = 6.67408e-11
        # Name of the body (string)
        self.id = id
        self.name = name
        # Initial position of the body (au)
        self.x0 = x0
        self.y0 = y0
        self.z0 = z0
        # Position (au). Set to initial value.
        self.x = self.x0
        self.y = self.y0
        self.z = self.z0
        # Initial velocity of the body (au/s)
        self.vx0 = vx0
        self.vy0 = vy0
        self.vz0 = vz0
        # Velocity (au/s). Set to initial value.
        self.vx = self.vx0
        self.vy = self.vy0
        self.vz = self.vz0
        # Mass of the body (kg)
        self.M = mass
        # Short name
        self.vector = vector

        self.totalforcex = totalforcex
        self.totalforcey = totalforcey
        self.totalforcez = totalforcez

# All Celestial Bodies

forcex = 0
forcey = 0
forcez = 0

Bodies = [
    CelBody(0, 'Sun', 1, 1, 1, 0, 0, 0, 1.989e30, 'sun', 0, 0, 0, 0, 0, 0),
    CelBody(1, 'Mercury', 1*au, 1, 1, 0, 29780, 0, 3.3e23, 'earth', 0, 0, 0, 0, 0, 0),
    ]

leftover_bin = []
templistx = []
templisty = []
templistz = []

for v in range(365242):
    for n in range(len(Bodies)):
        #Need to initialize the bodies

        planetinit = Bodies[n]

        for x in range(len(Bodies)):
            # Temporary lists and initial conditions
            planet = Bodies[x]

            if (planet == planetinit):
                pass

            else:
                rx = Bodies[x].x - Bodies[n].x
                ry = Bodies[x].y - Bodies[n].y
                rz = Bodies[x].z - Bodies[n].z

                r3 = (rx**2+ry**2+rz**2)**1.5
                gravconst = G*Bodies[n].M*Bodies[x].M
                fx = -gravconst*rx/r3
                fy = -gravconst*ry/r3
                fz = -gravconst*rz/r3


                # Make a temporary list of the total forces and then add them to get the resulting force
                templistx.append(fx)
                templisty.append(fy)
                templistz.append(fz)

        forcex = sum(templistx)
        forcey = sum(templisty)
        forcez = sum(templistz)
        templistx.clear()
        templisty.clear()
        templistz.clear()

        x = int(Bodies[n].x) + int(Bodies[n].vx) * dt
        y = int(Bodies[n].y) + int(Bodies[n].vx) * dt
        z = int(Bodies[n].z) + int(Bodies[n].vz) * dt

        Bodies[n].x = x
        Bodies[n].y = y
        Bodies[n].z = z

        vx = int(Bodies[n].vx) + forcex/int(Bodies[n].M)*dt
        vy = int(Bodies[n].vy) + forcey/int(Bodies[n].M)*dt
        vz = int(Bodies[n].vz) + forcez/int(Bodies[n].M)*dt

        Bodies[n].vx = vx
        Bodies[n].vy = vy
        Bodies[n].vz = vz

        t += dt




print(Bodies[0].name)
print(Bodies[0].x)
print(Bodies[0].y)
print(Bodies[0].z)


print(Bodies[1].name)
print(Bodies[1].x)
print(Bodies[1].y)
print(Bodies[1].z)

It should output something like the coordinates here, but then also a z coordinate: coordinate 1 (41.147123353981485, -2812171.2728945166) coordinate 2 (150013715707.77917, 2374319765.821534)

But it outputs the following:

Sun 0.0, 0.0, 0.0

Earth 149600000000.0, 0.0, 0.0

Note: The problem is probably in the for loops or in the rounding of the sum of the arrays but I am not sure.

Solution

picture - 1000 words

The direct errors in your code are

You compute the force in the wrong direction, it should be rx = b[n].x-b[x].x etc. or you need to remove the minus sign some lines later.
Your computation in single coordinates invites copy-paste errors as in
```
x = int(Bodies[n].x) + int(Bodies[n].vx) * dt
y = int(Bodies[n].y) + int(Bodies[n].vx) * dt
z = int(Bodies[n].z) + int(Bodies[n].vz) * dt
```
where in the y coordinate you still use vx. The intermediate rounding to integer values makes no sense, it only reduces the accuracy somewhat.

I changed your code to use numpy arrays as vectors, separated the acceleration computation from the Euler updates, removed the non-sensical rounding to integer values during the numerical simulation, removed unused variables and fields, removed intermediate variables for the force/acceleration computations to directly update the acceleration field, changed the loop to use the time to notice when a year (or 10) has passed (your code iterates over 100 years in 0.1day increments, was that intended?), ... and added Venus to the bodies and added code to produce images, result see above.

This spiraling is typical for the Euler method. You can easily improve that pattern by changing the Euler update to the symplectic Euler update, which means to update the velocity first and compute the position with the new velocity. This gives, with everything else the same, the image

day = 60*60*24
# Constants
G = 6.67408e-11
au = 1.496e11

class CelBody(object):
    # Constants of nature
    # Universal constant of gravitation
    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 (kg)
        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 Bodies

t = 0
dt = 0.1*day

Bodies = [
    CelBody(0, 'Sun', [0, 0, 0], [0, 0, 0], 1.989e30, 'yellow', 10),
    CelBody(1, 'Earth', [-1*au, 0, 0], [0, 29783, 0], 5.9742e24, 'blue', 3),
    CelBody(2, 'Venus', [0, 0.723 * au, 0], [ 35020, 0, 0], 4.8685e24, 'red', 2),
    ]

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

# loop over ten astronomical years
v = 0
while t < 10*365.242*day:
    # 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 = sum(rx**2)**1.5
            body.a += -G*other.M*rx/r3

    for n, planet in enumerate(Bodies):
        # use the symplectic Euler method for better conservation of the constants of motion
        planet.v += planet.a*dt
        planet.x += planet.v*dt
        paths[n].append( planet.x[:2].copy() )
        #print("%10s x:%53s v:%53s"%(planet.name,planet.x, planet.v))
    if t > v:
        print("t=%f"%t)
        for b in Bodies: print("%10s %s"%(b.name,b.x))
        v += 30.5*day
    t += dt

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()