Why does the moon spiral towards Earth in simulation?

Question

I am trying to simulate a somewhat realistic program where Earth and the moon can interact gravitationally with each other. Now the problem is that the moon keeps on spiraling towards Earth and I don't understand why.

This is my code:

from math import sin,cos,sqrt,atan2,pi
import pygame
pygame.init()

class Planet:
    dt = 1/100
    G = 6.67428e-11 #G constant
    scale = 1/(1409466.667) #1 m = 1/1409466.667 pixels
    def __init__(self,x=0,y=0,radius=0,color=(0,0,0),mass=0,vx=0,vy=0):
        self.x = x #x-coordinate pygame-window
        self.y = y #y-coordinate pygame-window
        self.radius = radius
        self.color = color
        self.mass = mass
        self.vx = vx #velocity in the x axis
        self.vy = vy #velocity in the y axis
        
    def draw(self,screen):
        pygame.draw.circle(screen, self.color, (self.x, self.y), self.radius)
    
    def orbit(self,trace):
        pygame.draw.rect(trace, self.color, (self.x, self.y, 2, 2))
        
    def update_vel(self,Fnx,Fny):
        ax = Fnx/self.mass #Calculates acceleration in x- and y-axis for body 1.
        ay = Fny/self.mass
        self.vx -= ((ax * Planet.dt)/Planet.scale)
        self.vy -= ((ay * Planet.dt)/Planet.scale)
        self.update_pos()
     
    def update_pos(self):
        self.x += ((self.vx * Planet.dt)) #changes position considering each body's velocity.
        self.y += ((self.vy * Planet.dt))
        
    def move(self,body):
        dx = (self.x - body.x) #Calculates difference in x- and y-axis between the bodies
        dy = (self.y - body.y)
        r = (sqrt((dy**2)+(dx**2))) #Calculates the distance between the bodies
        angle = atan2(dy, dx) #Calculates the angle between the bodies with atan2!
        if r < self.radius: #Checks if the distance between the bodies is less than the radius of the bodies. Uses then Gauss gravitational law to calculate force.
            F = 4/3 * pi * r
            Fx = cos(angle) * F
            Fy = sin(angle) * F
        else:  
            F = (Planet.G*self.mass*body.mass)/((r/Planet.scale)**2) #Newtons gravitational formula.
            Fx = cos(angle) * F
            Fy = sin(angle) * F
        return Fx,Fy

def motion():
    for i in range(0,len(bodies)):
        Fnx = 0 #net force
        Fny = 0
        for j in range(0,len(bodies)):
            if bodies[i] != bodies[j]:
                Fnx += (bodies[i].move(bodies[j]))[0]
                Fny += (bodies[i].move(bodies[j]))[1]
            elif bodies[i] == bodies[j]:
                continue
        bodies[i].update_vel(Fnx,Fny)
        bodies[i].draw(screen)
        bodies[i].orbit(trace)
        Fnx,Fny=0,0 

screen = pygame.display.set_mode([900,650]) #width - height
trace = pygame.Surface((900, 650))
pygame.display.set_caption("Moon simulation")
FPS = 150 #how quickly/frames per second our game should update.

earth = Planet(450,325,30,(0,0,255),5.97219*10**(24),-24.947719394204714/2) #450= xpos,325=ypos,30=radius
luna = Planet(450,(575/11),10,(128,128,128),7.349*10**(22),1023)
bodies = [earth,luna]

running = True
clock = pygame.time.Clock()

while running: #if user clicks close window
    clock.tick(FPS)    
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False
            
    screen.fill((0,0,0))
    pygame.Surface.blit(screen, trace, (0, 0))
    motion()

    pygame.display.flip() 

pygame.quit()

Once I get the earth-moon system to work I would like to expand on it and try having three bodies (the reason why there is so much of what otherwise would be "unnecessary" code)

I am open to suggestions or/and advices! Thanks

Short answer: rounding errors. Long answer: https://towardsdatascience.com/modelling-the-three-body-problem-in-classical-mechanics-using-python-9dc270ad7767 — MatBailie, Mar 08 '23 at 22:17
Cool project! My hunch is: truncation error, due to the chosen approximation. Rounding errors would not be my first guess, because these tend to be random/smeary, but usually not biased. (I only read the headline and had a glance at your code so far.) What method are you using to approximate the PDE? Euler-forward, Euler-backward, or some higher-order method? If you use a low-order method (like Euler), then the answer is definitly: truncation error. — Christian Fuchs, Mar 08 '23 at 22:22
@user2357112 What do you mean? When updating my velocities, I multiply dt with the acceleration — Hale, Mar 08 '23 at 22:37
@Hale If with semi-implicit you mean that you do a forward-update of the position (+) and a then a backward update of the velocity (-), then the spiral is probably due to the bias of your method. — Christian Fuchs, Mar 08 '23 at 22:53
@ChristianFuchs To be honest I am not entirely sure if this is the correct name of my method. I suppose this is the method I mentioned before https://en.wikipedia.org/wiki/Semi-implicit_Euler_method. All I know is that I begin with calculating the net force in the x- and y-axis separately on a body in order to get the acceleration so I can use it to calculate the velocity and distance. I do this with every body and update their positions accordingly. — Hale, Mar 08 '23 at 23:05
Note that if you specifically want to model orbits, then you can use something called a **symplectic integrator.** A symplectic integrator is designed to preserve the energy in the system, and it’s a common choice for celestial mechanics. — Dietrich Epp, Mar 12 '23 at 02:43

Christian Fuchs · Accepted Answer · 2023-03-12T03:00:32.057

You need to compute the forces in one go, not body-move by body-move.

Instead of computing the forces during the update loop, while shifting positions:

def motion():
    for i in range(0,len(bodies)):
        Fnx = 0 #net force
        Fny = 0
        for j in range(0,len(bodies)):
            if bodies[i] != bodies[j]:
                Fnx += (bodies[i].move(bodies[j]))[0]
                Fny += (bodies[i].move(bodies[j]))[1]
            elif bodies[i] == bodies[j]:
                continue
        bodies[i].update_vel(Fnx,Fny)
        bodies[i].draw(screen)
        bodies[i].orbit(trace)
        Fnx,Fny=0,0

settle the forces upfront:

def motion():
    force = [ (
        sum([(bodies[i].move(bodies[j]))[0] for j in range(0, len(bodies)) if i != j ]),
        sum([(bodies[i].move(bodies[j]))[1] for j in range(0, len(bodies)) if i != j ])
    ) for i in range(0,len(bodies)) ]
    for i in range(0,len(bodies)):
        Fnx = force[i][0]
        Fny = force[i][1]
        bodies[i].update_vel(Fnx,Fny)
        bodies[i].draw(screen)
        bodies[i].orbit(trace)
        Fnx,Fny=0,0

(I don't write in Python usually, so the style is not perfect.)

The following text is from a previous answer. It may be helpful, but it is not required to solve the problem asked; you may stop reading here.

You can reduce numeric truncation errors further with more elaborated methods, like Runge-Kutta. For that:

don't do update_vel and update_pos one-after the other, instead, try to write an update_state method which combines both simultaneously; important is that the left hand side of the equations is either the delta or the new state, the right hand side of the equations is the old state, exclusively (higher-order Runge-Kutta will have some intermediate states, fractions of Planet.dt)

If Runge-Kutta is too heavy to start with, consider MacCormack or Lax-Wendroff.

For a Lax-Wendroff'ish way, instead of:

    def update_vel(self,Fnx,Fny):
        ax = Fnx/self.mass
        ay = Fny/self.mass
        self.vx -= ((ax * Planet.dt)/Planet.scale)
        self.vy -= ((ay * Planet.dt)/Planet.scale)
        self.update_pos()

    def update_pos(self):
        self.x += ((self.vx * Planet.dt))
        self.y += ((self.vy * Planet.dt))

try this:

    def update_state(self,Fnx,Fny):
        ax = Fnx/self.mass
        ay = Fny/self.mass
        self.x += (self.vx * Planet.dt) - (ax/Planet.scale) * Planet.dt**2
        self.y += (self.vy * Planet.dt) - (ay/Planet.scale) * Planet.dt**2
        self.vx -= (ax/Planet.scale) * Planet.dt
        self.vy -= (ay/Planet.scale) * Planet.dt

@Hale after you wrote that you use some semi-implict Euler method, I have rewritten my answer; hope it is still useful! — Christian Fuchs, Mar 08 '23 at 23:00
Thank you so much for your answer and advice!!! As for the Runge-Kutta method, is there any other (perhaps easier) method to implement? It's quite complicated due to the fact that I haven't begun learning about differential equations yet nor how to solve them. — Hale, Mar 08 '23 at 23:20
Methods between Runge-Kutta and 1st-order Euler are e.g. MacCormack or Lax-Wendroff. Good luck! — Christian Fuchs, Mar 08 '23 at 23:34
when exchanging my code with above I get an orbit looking like an outward spiral instead of a circle shaped one. — Hale, Mar 09 '23 at 10:28
@Hale you're right; there were two problems, one in my answer (I fixed that, only `**2` instead of `**2/2` for the 2nd order term), and another one in your `motion` method; now I compute all the forces before the loop, not during the loop, that resolves the problem. — Christian Fuchs, Mar 10 '23 at 00:37
@Hale I realized that I get a circular path by fixing the `motion` method only, regardless of the `update_X` method (at least for short simulations), which again shows that your equations are not wrong themselves. So I have again updated my answer for that, with emphasis on the `motion` method — Christian Fuchs, Mar 10 '23 at 00:47
It works!!! Now it is a circular path with the orbit looking as it should! Thank you so so much for your time and help!!! I deeply appreciate it! — Hale, Mar 11 '23 at 13:45

Why does the moon spiral towards Earth in simulation?

1 Answers1