鐘擺模擬

Question

我想在Pygame中編寫簡單的鐘擺模擬。重點是我試圖直接模擬擺錘上的力（重力和張力），而不是求解描述運動的微分方程。首先我寫了一個函數，它獲得一個矢量，以某個角度旋轉軸系，並將這個矢量的分量返回到新的旋轉軸系統;這個函數的代碼很好，它按預期工作。

模擬中的每個勾號I將重力矢量旋轉擺錘和繩索之間的角度，並獲得新組件 - 一個在繩索方向，一個與其正交。繩索方向上的張力和分量相互抵消，因此只有正交分量是重要的。計算完後，我將加速度矢量旋轉回正常座標系並進行積分。但是，由此產生的行爲並非如預期的那樣。可能是什麼原因？

這是代碼：

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。