Pendulum simulation

Question

I am trying to write simple pendulum simulation in Pygame. The point is that I am trying to simulate directly the forces on the pendulum (gravity and tension) rather than solving the differential equation that describes the motion. First I wrote a function that get a vector, rotate the axis-system by some angle, and return this vector's components in the new, rotated axis-system; the code of this function is fine and it works as expected.

Each tick of the simulation I rotate the gravity vector by the angle between the pendulum and the rope, and get the new components - one in the direction of the rope, and one is orthogonal to it. the tension the and component in the direction of the rope cancelling each other, so only the orthogonal component is important. After I calculate it, I rotate the acceleration vector back to the normal coordinates system, and integrate. However, the resulting behavior is not as intended. What can be the reason?

This is the code:

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

Thank you.

Answer 1

There are a bunch of problems here. I'll fix a few and leave a few for you.

What I fixed was: 1) since you've imported numpy, you should use it, and write things in terms of the vectors; 2) it's an unreasonable demand on yourself to write everything and have it work immediately; so you need to plot intermediate results, etc, like here I plot a as well, so you can see whether it makes sense; 3) your whole "rotation" approach is confusing; instead think of component parts; which I calculate here directly (it's shorter, easier to read and understand, etc); 4) in all simulations where you use a time step, you should explicitly use dt so you can change the timestep without changing other parameters.

Now if you watch it you can see it looks almost reasonable. Notice though that the acceleration never goes upward, so the ball just falls while it oscillates. The reason for this is that you did not include the tension of the rope into the forces on the ball. I'll leave that part to you.

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

If you want to do a simulation, I think you're doing it the hard way. I'd start with the equation of motion, see Equation 20, here. The dots mean differentiate with respect to time---so the equation is d^2/dt^2 \theta = ... Then you should implement a finite differences scheme in the time direction, and step through time. At each step (label with i), you can calculate the x and y coordinates of the bob, based on the length of the pendulum and \theta_i. Check out the wiki article on finite differences.