One of the things that really made physics “click” for me was learning how to numerically simulate the systems I learned about in class. It’s one thing to write down the equations of motion for a certain problem – it’s another to actually visualize how they look like on a computer. It’s the missing link between the abstract and the real.

In this post, I’ll solve the problem of the simple pendulum – a ubiquitous toy problem in physics – and show how to simulate it on a computer. To understand the concepts in this post, you should have a basic knowledge of calculus and be familiar with programming logic.

The first step to simulating any sort of physical problem is to write out the analytic equations of motion. We will use the Lagrangian method, a particularly powerful way to analyze complex systems. Using it for the case of a simple pendulum is probably overkill, but I like the elegance in its approach. Our first step is to define our problem. A schematic of a simple pendulum is presented below:

The problem statement is then this: **Find the equation of motion for a pendulum of length l and mass m. **The Lagrangian approach to doing this is as follows:

- Identify the generalized coordinate of the system, .
- Write the Lagrangian, where is the kinetic energy and is the potential energy, in terms of the generalized coordinates .
- The equations of motion for the system will then be given by the Euler-Lagrange equation, .

There will always be as many generalized coordinates as there are degrees of freedom in your system. The way I like to think about generalized coordinates is to imagine I’m on the phone with someone and I need to describe the state of my system to them. What is the minimum set of information I can tell them so that they can completely replicate my system? In the case of a simple pendulum, I could tell them the and positions of the mass, of course, but I could also just tell them the angle the pendulum makes with the vertical. One is fewer than two, and so the generalized coordinate for our simple pendulum is simply , corresponding to one degree of freedom.

The next step is to write out our energies in terms of the angle . The kinetic energy is easily written by remembering the relation for circular motion and noting that in our case, the radius is simply the length of the pendulum:

For the potential energy, we define our “zero” as when the mass is perfectly vertical, i.e. when . Then the “height” above the zero will just be , so the gravitational potential energy is:

We can now write the Lagrangian:

and apply the Euler-Lagrange equation to get our equation of motion:

Now we have our final analytic equation of motion. Note that the mass terms cancel out, suggesting that the motion of a pendulum is independent of its mass. At this stage, many introductory physics courses will take the small-angle approximation in order to obtain the equation for simple harmonic motion, which can be solved analytically. Without this approximation, there is no analytic solution for the simple pendulum, but that won’t bother us here, since we are seeking to solve it numerically. We have a second-order ordinary differential equation, which we can write as two first-order ordinary differential equations:

To numerically solve a system, the first step is discretizing it. The usual way to do this is by writing out the Taylor series for a continuous function and truncating it at some term. In our case, we want to find an approximation to given our known initial condition . Expanding about , gives, to the first order:

Rewriting, we have:

and plugging this into our system of differential equations, we get:

There is a question though – for the right hand side, at what time value should we compute and ? This is a non-trivial thing to decide, but for now we will compute these values at the past value , else we end up with multiple unknowns to solve for. Moving all our “knowns” to the right-hand side, we have:

We now have a method for determining and given the initial conditions and where . Using and for and and and for and , we can translate this to programming logic as follows:

omega = omega_old - (g/l)*sin(theta_old)*tau theta = theta_old + omega_old*tau

where we have defined to be a constant number. This is our update equation, which we will loop through every time we want to “step forward” in time to obtain a future value for . The only things we need to specify ahead of time, aside from the constants, are the initial conditions. The relevant code, written in Python, looks like this:

import numpy as np from math import * def simplePendulumSimulation(theta0,omega0,tau,m,g,l, numSteps,plotFlag): # initialize vectors time_vec = [0]*numSteps theta_vec = [0]*numSteps omega_vec = [0]*numSteps # set initial conditions theta = theta0 omega = omega0 time = 0 # begin time-stepping for i in range(0,numSteps): omega_old = omega theta_old = theta # update the values omega = omega_old - (g/l)*sin(theta_old)*tau theta = theta_old + omega_old*tau # record the values time_vec[i] = tau*i theta_vec[i] = theta omega_vec[i] = omega

Running this code for initial conditions and , we immediately see a problem, though:

The total energy is increasing without bound! In other words, the energy of the system is not conserved. You can see this by the ever-increasing displacement of our pendulum, which paradoxically swings higher and higher with every rotation. The problem lies in our crude first-order approximation we made earlier when discretizing our solution. After all, our equation of motion is nonlinear, making linear approximations especially ineffective at solving it. In lieu of using higher-order methods, we can just apply a simple modification called the Euler-Cromer method, which is guaranteed to conserve energy. We simply use the updated value of the angular velocity as it becomes available:

omega = omega_old - (g/l)*sin(theta_old)*tau theta = theta_old + omega*tau

You can see that in the second line has been simply changed to . Running the simulation again, we get the results:

The pendulum now properly returns to its initial displacement after each cycle. We see that although the total energy still fluctuates a bit (due to our first-order approximation), it is now properly bounded. The Python script to perform the simulation using the Euler-Cromer method and generate the images is given here: pendulum.py