摆振模拟

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

如果你想做一个模拟，我觉得你可能走了弯路。首先，我建议你从运动方程开始，具体可以参考第20个方程。这里的点表示对时间进行微分——所以这个方程是d^2/dt^2 \theta = ...。接下来，你应该在时间上实现一个有限差分的方法，然后一步一步地推进时间。在每一步（用i标记），你可以根据摆的长度和\theta_i来计算摆锤的x和y坐标。你可以查看一下维基百科关于有限差分法的文章。

Answer 2

这里有很多问题。我会解决一些，留下一些给你自己去处理。

我解决的内容包括：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()