Simulation of Pendulum: VPython Tutorial 3 (Visual Python)

Simulation of Pendulum
Simulation of Pendulum using VPython. (Credit: Tech for Curious)

A pendulum is the first instrument which comes to mind when someone thinks about Physics. The pendulum is a simple instrument which consists of a bob attached to a string whose other end is tied to a pivot at some height. When the pendulum is left on itself, the bob hangs down with stretched string and remain stationary. If someone gives bob slight push, the pendulum starts moving and swings back and forth. The push provides the pendulum a kinetic energy which helps in the motion of the pendulum.

One can set the pendulum in motion by raising the height of pendulum also by moving it sideway while the string is still in stretched position. In doing so, one actually raises the height of bob which supplies the pendulum a gravitational potential energy. This back and forth motion of the pendulum is called oscillation. Thus, you can oscillate a pendulum by either giving is a push or raising the bob to some height.

When the bob reaches the maximum height, the velocity becomes zero. At this point, the kinetic energy of pendulum is zero and potential energy is maximum. When the bob reaches the lowest height, then the velocity is maximum. At this point, the kinetic energy is maximum and potential energy is minimum. The total energy of pendulum is the sum of potential and kinetic energy.

An object of mass m moving with velocity v has a kinetic energy given as

(1)   \begin{equation*} K.E.=\frac{1}{2}mv^2 \end{equation*}

Whereas, the potential energy gained by an object when it is moved upward by a height h is given as

(2)   \begin{equation*} P.E.=mgh \end{equation*}

where g is acceleration due to gravity (9.8 m/s2).

Diagram of Pendulum
The schematic Diagram of pendulum showing gravitational force and its components. (Credit: Tech for Curious)

The force which helps the motion of the pendulum is the gravitation force which is equal to mg and acts downward on the bob. Suppose at any particular time, the angle between the string and the vertical direction is θ. The force mg can be resolved into two components mg(cosθ) and mg(sinθ). The component mg(cosθ) is balanced by the tension of the string T and the component mg(sinθ) is responsible for the motion of the pendulum.

(3)   \begin{equation*} T=mg\cos{\theta} \end{equation*}

(4)   \begin{equation*} ma=-mg\sin{\theta} \end{equation*}

(5)   \begin{equation*} l\frac{d^2\theta}{dt^2}=-g\sin{\theta} \end{equation*}

(6)   \begin{equation*} \frac{d^2\theta}{dt^2}=-\frac{g}{l}\sin{\theta} \end{equation*}

where l is the length of the pendulum. We can find the angle θ, using the numerical method as we did earlier in last two simulations for finding positions. Let us represent angular velocity and acceleration as

(7)   \begin{equation*} \omega=\frac{d\theta}{dt} \end{equation*}

(8)   \begin{equation*} \alpha=\frac{d^2\theta}{dt^2}=-\frac{g}{l}\sin{\theta} \end{equation*}

Suppose, at time t, the pendulum is at some angle θ(t) and moving with angular velocity ω(t),  then the angular acceleration at that time t is

(9)   \begin{equation*} \alpha(t)=-\frac{g}{l}\sin(\theta(t)) \end{equation*}

The angular position after time interval dt can be found as

(10)   \begin{equation*} \theta(t+dt)=\theta(t)+\omega(t)dt \end{equation*}

We can find the angular velocity after time interval dt as

(11)   \begin{equation*} \omega(t+dt)=\omega(t)+\alpha(t)dt \end{equation*}

Thus, we can see that just like the linear position, we can calculate the angular position also but we would like to find the position of the pendulum in terms of x,  y and z coordinates. Usually, we represent the motion of the pendulum in x-y plane and the z coordinate of the pendulum remains constant throughout the motion. The directions towards the right is +ve and left is -ve whereas in vertical directions, the upward direction is +ve and downward is -ve. According to the diagram, the coordinates of bob are equal to

(12)   \begin{equation*} x=HP=l\sin{\theta} \end{equation*}

(13)   \begin{equation*} y=Pivot.y-l\cos{\theta} \end{equation*}

The summary of the above concept is that

  • Calculate α(t) from θ(t)
  • Update θ(t)
  • Update ω(t)
  • Calculate x(t) and y(t).

We will be writing code in VPython for our simulation, if you are doing VPython first time you can first go through the Simulation of Bouncing Ball or the Simulation of Mass-Spring System although it is not prerequisite. We need to open VIDLE to write the code although you can write code in simple text editor also but it is more convenient in VIDLE.

We will first import visual python module and define our display.

from visual import*
display(width=600,height=600,center=vector(0,12,0),background=color.white)

This code will set the display size to 600 x 600, center at y = 12 and x = z =0 and the color of the background will be white. The horizontal direction is x-axis, the vertical direction is y-axis and the direction perpendicular to the screen is z axis. You can run the code by pressing F5 or clicking Run>Run Module, you will see only a white window of size 600 x 600. Now let us create pendulum.

from visual import*
display(width=600,height=600,center=vector(0,12,0),background=color.white)
g=9.8 # acceleration due to gravity
bob=sphere(pos=vector(5,2,0),radius=0.5,color=color.blue)
pivot=vector(0,20,0)
roof=box(pos=pivot,size=vector(10,0.5,10),color=color.green)
rod=cylinder(pos=pivot,axis=bob.pos-pivot,radius=0.1,color=color.red)

Run this code

Edit this code

The third line defines the acceleration due to gravity. The line after # is a comment which is ignored by the compiler. The comment helps in understanding the code later. The fourth line defines the bob of the pendulum which is a sphere of radius 0.5 unit and the initial position of the bob is (5,2,0). The pivot is a point where the string of pendulum is attached to the roof or rigid support, we need this to define the length of string with respect to the position of bob. Next, the roof is just a box whose central position is present at the position of the pivot. The length and breadth of the roof are 10 units each whereas the thickness is 0.5 unit in the y direction.

Stationary Pendulum
The objects like bob, string, and roof are created with some initial position. (Credit: Tech for Curious)

The string or rod of the pendulum is defined as a cylinder with radius 0.1 unit. One end of the cylinder is present at the position of the pivot and another end is present at the position of bob. The other end can be described by the axis which is not the actual position of the other end instead it is the relative position of the other end with respect to the first end. Thus, we give the difference of the two end’s positions as the value to the axis of the cylinder. After running the code, we see that the pendulum is created which is attached to a roof but the bob is stationary.

Now, we would like to update the position of bob according to the equations described above.

from visual import*
display(width=600,height=600,center=vector(0,12,0),background=color.white)
g=9.8 # acceleration due to gravity
bob=sphere(pos=vector(5,2,0),radius=0.5,color=color.blue)
pivot=vector(0,20,0)
roof=box(pos=pivot,size=vector(10,0.5,10),color=color.green)
rod=cylinder(pos=pivot,axis=bob.pos-pivot,radius=0.1,color=color.red)
t=0 # time 
dt=0.01 # time interval 
l=mag(bob.pos-pivot) # length of pendulum
cs=(pivot.y-bob.pos.y)/l # calculation of cos(theta) 
theta=acos(cs) # angle with vertical direction
vel=0.0 # angular velocity
while (t<100):
 rate(100) # maximum 100 calculations per second
 acc=-g/l*sin(theta) # updating of angular acceleration
 theta=theta+vel*dt # updating of angular position
 vel=vel+acc*dt # updating of angular velocity
 bob.pos=vector(l*sin(theta),pivot.y-l*cos(theta),0) # cal. position
 rod.axis=bob.pos-rod.pos # updating other end of rod of pendulum
 t=t+dt # updating time

Run this code

Edit this code

We have defined many new variables like t, dt and l which specify the time, time interval and length of pendulum respectively. The mag() function determines the magnitude of any vector which we have used to determine the distance between the position of the bob and the pivot. This actually gives the length of the string of the pendulum. We have calculated the cosθ from the formula

(14)   \begin{equation*} \cos{\theta}=\frac{Pivot.y-y}{l} \end{equation*}

and by taking the cos inverse, we can find the angle θ at the starting point. We have calculated θ because our equation required angular position θ not the position (x,y). The initial angular velocity is taken as zero.

Now, the initial angular position and angular velocity are defined. We can update the angular position and velocity inside a loop to get a moving pendulum. The codes written inside the loop are executed till t is less than 100. We have first calculated the angular acceleration from the angular position theta using the equation described above. Since the angular position is not required to calculate anything else so we will update the angular position using the angular velocity. After this, the angular velocity is updated using angular acceleration. Note that we should not update angular velocity before updating angular position because in that case, we will be using velocity at time (t+dt) to update position at the time (t) which give rise to an error in prediction of position. We should use velocity at time t to update position at time t.

Simulation of Pendulum
The pendulum is moving back and forth in oscillatory motion which is also called harmonic motion. (Credit: Tech for Curious)

Once we get the updated angular position θ, we can calculate the actual position in terms of x and y coordinates. Now, the position of the bob is updated but the position of another end of the string is still same so, we should update the axis of string also. In this way, the other end of the string will always be moving along the bob. In the last, we updated the time by adding small time interval dt.

After running the code you will see that the pendulum is moving as we expected. The pendulum starts from the right-hand side at some height and as it goes down, it gains velocity. The pendulum gets maximum velocity when it is at the bottom position which is also called equilibrium position. Since the bob still has kinetic energy and inertia, the bob start climbing in the left direction till it reaches the maximum height which should be same as the height it started with. Then the pendulum reverses back and retraces its path.

You can also set the pendulum in motion by starting at the bottom most position and giving it some initial velocity. I have made a video of the above tutorial in which I have shown all steps and shown that the pendulum can be set in motion by giving initial velocity also. If you have any query or suggestion, please fill free to use the comment box. Your suggestions help us to decide future tutorials. If you like this tutorial, please share with someone who is interested in visualizing physics. You can follow us on facebook and twitter. You can subscribe us for Email Notification also to get an email whenever we publish a new post.


4 thoughts on “Simulation of Pendulum: VPython Tutorial 3 (Visual Python)”

  1. Greetings,

    Thank you for sharing this tutorial. As I ran the program, I decided to graph the motion of the pendulum. I noticed an error but I don’t know the solution. Maybe if I address the error, someone will know how to fix it.

    Error:
    If you run (time vs. theta), you will notice the amplitude of the wave increasing as time increases.

    If you graph (time vs. bob.pos.y) when the initial velocity is 0.0, then again, I got the same result. The bob is dropped from an initial height of 2.0 and yet the bob climbs above this value for each period of the bob. This doesn’t follow the conservation of energy.

    The code you have simulates a form of “negative gravity” yet I’m not sure how.

    1. Hello Adam,
      Thank you for your comment. You are right, if you plot the total energy of the pendulum, it may seem to vary with time. The error lies in the time-step, we are using, theoretically, the time-step dt should be infinitesimal. But we are using dt as 0.01 which is a finite value. If you make dt smaller, you will see that the error which you are observing will become less.

      This problem of energy conservation, I have discussed in the next post where I have plotted the phase diagram of motion. There we can see that phase diagram do not retrace its path even for periodic motion but for smaller value, to dt it starts retracing its path more closely.

      I hope this answer satisfies your doubt. Please keep visiting and giving your valuable comments.

  2. Hello,

    I was trying to model a swinging motion for an assignment, so I tried changing the values of this code. But somehow, when I changed the initial position of the bob, everything just disappeared and I am not sure why. When I reset the initial position of bob back to (5,2,0) like you had in the video everything went back to normal. Can you please help explain why this happens?

    1. Hello Karina Patricia,
      First of all, thank you for your comment and pointing out the mistake. I have checked the code and found that there is a mistake in the convention of coordinates. The y coordinate of bob should be defined as Pivot.y-lcos(theta) rather than just HE because E is not located at origin and y of bob should be measured from the origin. Similarly cos(theta) = (Pivot.y-bob.y)/l. Now, the simulation is running fine.

Leave a Reply

Your email address will not be published. Required fields are marked *