摆锤模拟

Question

我正在尝试在 Pygame 中编写简单的钟摆模拟。关键是我试图直接模拟摆锤上的力（重力和张力），而不是求解描述运动的微分方程。首先，我编写了一个函数，它获取一个向量，将轴系旋转某个角度，然后在新的旋转轴系中返回该向量的分量；这个函数的代码很好，它按预期工作。

模拟的每个刻度，我将重力矢量旋转钟摆和绳索之间的角度，并获得新的组件 - 一个在绳索的方向上，一个与它正交。绳索方向上的张力和分量相互抵消，因此只有正交分量很重要。计算后，我将加速度矢量旋转回法线坐标系，并积分。但是，由此产生的行为并不符合预期。可能是什么原因？

这是代码：

from __future__ import division
import copy
import pygame
import random
import math
import numpy as np
import time

clock = pygame.time.Clock()
pygame.init()
size = (width, height) = (600,500)
screen = pygame.display.set_mode(size)

def rotate(vector,theta):
    #rotate the vector by theta radians around the x-axis
    Vx,Vy = vector[0],vector[1]
    cos,sin = math.cos(theta),math.sin(theta)
    newX,newY = Vx*cos-Vy*sin, Vy*cos+Vx*sin #the newX axis is the result of rotating x axis by theta
    return [newX,newY]

class pendulum:
    def __init__(self,x,y,x0,y0):
        self.x = x
        self.y = y
        self.x0 = x0
        self.y0 = y0
        self.velocity = [0,0]
        self.a= [0,0]
        self.angle = 0
    def CalcForce(self):
        self.angle = math.atan2(-(self.y-self.y0),self.x-self.x0)
        gravity = rotate(g,self.angle)
        self.a[1]=gravity[1]
        self.a[0] = 0 #This component is cancelled by the tension
        self.a = rotate(self.a,-self.angle)
    def move(self):
       #print pylab.dot(self.velocity,[self.x-self.x0,self.y-self.y0])
        self.velocity[0]+=self.a[0]
        self.velocity[1]+=self.a[1]
        self.x+=self.velocity[0]
        self.y+=self.velocity[1]
    def draw(self):
        pygame.draw.circle(screen, (0,0,0), (self.x0,self.y0), 5)
        pygame.draw.line(screen, (0,0,0), (self.x0,self.y0), (int(self.x), int(self.y)),3)
        pygame.draw.circle(screen, (0,0,255), (int(self.x),int(self.y)), 14,0)
g = [0,0.4]
p = pendulum(350,100,300,20)

while 1:
    screen.fill((255,255,255))
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
    p.CalcForce()
    p.move()
    p.draw()
    clock.tick(60)
    pygame.display.flip()

谢谢。

Answer 1

这里有很多问题。我会修几个，留几个给你。

我修复的是：1）既然你已经导入了 numpy，你应该使用它，并用向量来写东西； 2）对自己的要求是不合理的，把所有的东西都写出来并立即工作；所以你需要绘制中间结果等，比如这里我也绘制a，所以你可以看看它是否有意义； 3）您的整个“轮换”方法令人困惑；而是考虑组成部分；我在这里直接计算（它更短，更容易阅读和理解等）； 4) 在使用时间步长的所有模拟中，您应该明确使用 dt，这样您就可以在不更改其他参数的情况下更改时间步长。

现在，如果您观看它，您会发现它看起来几乎是合理的。请注意，尽管加速度永远不会向上，所以球只是在振动时下落。原因是您没有将绳索的张力包括在球上的力中。我会把这部分留给你。

import pygame
import math
import numpy as np

clock = pygame.time.Clock()
pygame.init()
size = (width, height) = (600,500)
screen = pygame.display.set_mode(size)


class pendulum:
    def __init__(self,x,y,x0,y0):
        self.x0 = np.array((x0, y0))
        self.x = np.array((x, y), dtype=float)
        self.v = np.zeros((2,), dtype=float)
        self.a = np.zeros((2,), dtype=float)
    def CalcForce(self):
        dx = self.x0 - self.x
        angle = math.atan2(-dx[0], dx[1])
        a = g[1]*math.sin(angle)  # tangential accelation due to gravity
        self.a[0] = at*math.cos(angle)
        self.a[1] = at*math.sin(angle)
    def move(self):
        #print np.dot(self.a, self.x-self.x0) #is a perp to string?
        self.x += dt*self.v
        self.v += dt*self.a
    def draw(self):
        pygame.draw.circle(screen, (0,0,0), self.x0, 5)
        pygame.draw.line(screen, (0,0,0), self.x0, self.x.astype(int),3)
        pygame.draw.circle(screen, (0,0,255), self.x.astype(int), 14,0)
        pygame.draw.line(screen, (255, 0, 0), (self.x+200*self.a).astype(int), self.x.astype(int), 4)
dt = .001
g = [0,0.4]
p = pendulum(350,100,300,20)

while 1:
    screen.fill((255,255,255))
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
    for i in range(100): # don't plot every timestep
        p.CalcForce()
        p.move()
        p.draw()
    clock.tick(60)
    pygame.display.flip()

Answer 2

如果你想做一个模拟，我认为你做的很艰难。我将从运动方程开始，请参阅Equation 20，在这里。点的意思是关于时间的微分---所以方程是d^2/dt^2 \theta = ... 然后你应该在时间方向上实现一个有限差分方案，并逐步通过时间。在每一步（标有i），您可以根据摆锤的长度和\theta_i 计算鲍勃的x 和y 坐标。查看wiki article 的有限差分。