Python pendulum

972 views Asked by At

I'm trying to deduce the pendulum equation out of the Verlet algorithm. This is what I've done with Python using r as the angle of the pendulum, w as the speed and a as the angular acceleration. h is the increase of time, and my main problem is that if I introduce an h lower than 2 the program just show me the same number again and again.

from math import sin, cos, pi
from pylab import plot,show
h=int(input("h: "))
q=int(input("Angle: "))
l=int(input("Pendulum lenght: "))
r=q*pi/180 #Change the angle from degrees to rad
w=0
g=9.81
a=(-g/l)*sin(r)
y=0
xx=[]
yy=[]
while y>=0:
    r= r+h*w+ (h**2)/2*a
    x= cos(r) #Cartesians
    y= sin(r)
    w= w+ ((h/2)*a)
    a=(-g/l)*sin(r)
    w= w+((h/2)*a)
    xx.append(x)
    yy.append(y)
plot (xx,yy)
show()
2

There are 2 answers

0
xnx On BEST ANSWER

For any physically-reasonable pendulum, your h should be less than 1 (seconds) to model this properly. But you're casting to int which rounds down, so you get h=0 (no time step). Instead, use float:

h=float(input("h: "))
q=float(input("Angle: "))
l=float(input("Pendulum length: "))
2
aaaantoine On

I don't know the first thing about physics -- xnx likely has your answer.

But, can I recommend that you take a functional approach to programming?

Let's start by encapsulating the formula for a:

from math import sin, cos, pi
from pylab import plot,show

def getA(pLen, r):
    g = 9.81
    return (-g / pLen) * sin(r)

Then we need to make the rest of your program a function so we can do something very useful in a moment.

def main():
    h=int(input("h: "))
    q=int(input("Angle: "))
    l=int(input("Pendulum lenght: "))
    r=q*pi/180 #Change the angle from degrees to rad
    w=0
    a=getA(l, r)
    y=0
    xx=[]
    yy=[]
    while y>=0:
        r= r+h*w+ (h**2)/2*a
        x= cos(r) #Cartesians
        y= sin(r)
        w= w+ ((h/2)*a)
        a=getA(l, r)
        w= w+((h/2)*a)
        xx.append(x)
        yy.append(y)
    plot (xx,yy)
    show()

# call your main function only when running directly
if __name__ == '__main__':
    main()

Now, assuming your file is named pendulum.py, jump to your console and open the python shell:

>python
Python 2.7.6 (default, Nov 10 2013, 19:24:18) [MSC v.1500 32 bit (Intel)] on win
32
Type "help", "copyright", "credits" or "license" for more information.
>>> import pendulum
>>> pendulum.getA(5, 50)
0.5147794629671083
>>>

You can now desk check your getA() function to make sure it's returning the right output for the given inputs.

Repeat for all other operations inside the program, and you'll be able to isolate your problem further, which is an important step in solving it.