模拟中为什么月球会朝地球螺旋而去？

Question

模拟中为什么月球会朝地球螺旋而去？

pythonpygamesimulationgame-physicsphysics

5

我正在尝试模拟一个相对真实的程序，其中地球和月球可以通过引力相互作用。现在问题是月球不断向地球螺旋，并且我不明白为什么。

这是我的代码：

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

一旦我让地月系统运作起来，我想要扩展它并尝试三体问题（这就是为什么有那么多本来会是“不必要”的代码的原因）

我愿意听取建议和/或意见！谢谢

- Hale

2

简短回答：舍入误差。长篇回答请参考：https://towardsdatascience.com/modelling-the-three-body-problem-in-classical-mechanics-using-python-9dc270ad7767 - MatBailie

1

很酷的项目！我的直觉是：由于所选的近似值，导致截断误差。舍入误差不会是我的第一选择，因为这些误差往往是随机/模糊的，但通常不是有偏的。（我只看了标题并浏览了你的代码）你用什么方法来近似PDE？欧拉前向、欧拉后向还是某种高阶方法？如果你使用低阶方法（如欧拉），那么答案肯定是：截断误差。 - Christian Fuchs

哦，你懂。我不知怎么把代码读错了。 - user2357112

@Hale，这是一个兔子洞，无论如何，祝你探索愉快！ - Christian Fuchs

1

请注意，如果您想特别建模轨道，则可以使用称为辛积分器的工具。辛积分器旨在保持系统中的能量，并且是天体力学的常见选择。 - Dietrich Epp

显示剩余4条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Christian Fuchs · Accepted Answer

你需要一次性计算出所有力的作用，而不是逐个物体地进行计算。

与其在更新循环中计算力量，同时移动位置：

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

提前解决力量问题：

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

（我通常不写Python，所以风格可能不太完美。）

以下文字来自之前的回答。它可能有帮助，但不是解决问题所必需的；您可以在此停止阅读。

您可以使用更详细的方法（如龙格-库塔）进一步减少数字截断误差。为此：

不要依次执行update_vel和update_pos，而是尝试编写一个update_state方法，同时将两者结合起来；重要的是方程式的左边是增量或新状态，方程式的右边是旧状态，仅限于这种情况（高阶龙格-库塔会有一些中间状态，即Planet.dt的分数）

如果龙格-库塔对您来说太重了，请考虑使用MacCormack或Lax-Wendroff。

对于类似于Lax-Wendroff的方法，而不是：

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

试一下这个：

    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