One of the more interesting things you can do with a computer is to write a program based around a simple physical law, such as f = m × a (force equals mass times acceleration), and use that program to see the complex and beautiful behavior of a physical system. For example, you could write your own simulator for celestial bodies, and use that simulator to explore the motion of the solar system, as well as of the earth and moon. This lecture will go over the Newtonian physics you’d need, and give techniques for simulating those physics with a computer program.
I would like to write a program to make a ball bounce on the screen under the influence of gravity. It’s a complicated problem: I have to take into account acceleration due to gravity, keeping track of positions and velocities, a floor and walls that the ball will bounce off of, and some graphics to draw as well. How fast should the simulation be, what are the units to use, and do I know the gravitational constant measured in units of pixels/sec2 (pixels per second-squared)?
Having identified all those issues, I admit to myself that I’m just not clever enough to solve all those subproblems at once.Here’s something I can handle. I’ll start with just making the ball move up and down at a constant speed.
Before we get to the program, a note about physics definitions. Velocity and speed are not the same. Speed says how fast an object is moving, but it says nothing about the direction in which it’s moving. Velocity combines speed with a direction. So “55 miles per hour” is a speed, and “heading northeast at 55 miles per hour” is a velocity. If you’ve ever seen vectors, you’ll recognize velocity as a vector and speed as the magnitude of the velocity.
Our program will keep track of the velocity, which we decompose into horizontal (x) and vertical (y) components. The sign of each component (+ or –) will tell us the direction in each component.
At this point in the course, you’ve done a couple of animations, so you can probably figure out how this program works. Because you have already written the pong game in Lab Assignment 1, you should be able to understand quite easily the code to make the ball bounce off the floor and ceiling.
I’d like to make the animation a bit smoother. In the last example, we drew frames at a rate of 20 frames per second. But my computer is fast enough to easily draw 60 frames per second.
The problem with increasing the frame rate is that each time draw_frame is called, the ball moves a constant amount:
y = y + v_ySo if draw_frame is called 60 times per second, instead of 20, the ball will move three times as fast. (This was a problem in early PC computer games: the game would draw as many frames as possible per second on the computer the game was written on. Then a faster computer would come out, and the evil ghosts would run far too fast. Many computers started to have a “turbo boost” button on the front that when pressed, let the computer run at full speed.)
The solution is to compute how much time has passed since the last time draw_frame was called, and multiply that amount of time by the velocity to find out how much to change the position of the ball. We call this amount the timestep. The simplest way to get an approximate value for the timestep is to divide 1 by the frame rate. (Making sure, of course, that we do floating division, not integer division.)
While we’re at it, let’s refactor the code and move some of the drawing and computation code into their own functions, simplifying the body of the main function.
If you look at the body of the compute_next_position function, it computes the new position as position + velocity * timestep. Let’s verify that the units work out. Position is in pixels, velocity is in pixels per second, and the timestep is in seconds. The units for the expression work out to
$\mbox{pixels} + \left(\displaystyle\frac{\mbox{pixels}}{\mbox{second}} \times \mbox{seconds}\right)$ ,
which is what we want.
In the original program, the ball’s speed was 4 pixels per frame, so at 20 frames per second, the ball’s speed was 80 pixels per second. That’s why, now that the velocities are in terms of pixels per second in the new version, I initialized v_y to —80.
Now we’ve solved the simpler problem of bouncing the ball with constant speed adequately, and we’re ready to tackle the original problem, where the ball accelerates due to gravity. We’ll remove the ceiling, letting gravity bring the ball back to the floor.
What is acceleration? It’s how the velocity changes over time. Think of it this way. The velocity measures how the ball’s position changes over time, and the acceleration measures how the ball’s velocity changes over time.
Suppose that the ball is initially at y = 0 meters and its velocity is 1 meter per second. Let’s use a timestep of 1 second. So, after 1 second the ball is at y = 1 meter, after 2 seconds it’s at y = 2 meters, after 3 seconds it’s at y = 3 meters, and so on. But now let’s accelerate the ball at 1 meter per second squared. That means in each second, the velocity will increase by 1 meter per second, so that the increase in velocity is 1 meter per second, per second; or 1 meter per second squared. We start with the ball once again at y = 0 meters and a velocity of 1 meter per second. After 1 second, the ball is at y = 1 meter, but now its velocity is 2 meters per second. So after 2 seconds, the ball is at y = 3 meters, having moved an additional 2 meters in that second, and its velocity is now 3 meters per second. After 3 seconds, the ball is at y = 6 meters, having moved an additional 3 meters, and its velocity is now 4 meters per second. After 4 seconds, the ball is at y = 10 meters with a velocity of 5 meters per second. And so on. As the ball accelerates, the change in its position from the beginning to the end of each timestep increases.
Now we see that we have to address two issues:
We can solve the second issue by changing all coordinates in the program to be measured in meters. Then we can have a constant PIXELS_PER_METER that we use to actually draw the ball and floor in the correct location. Let’s choose a value of 10.0 pixels per meter, and we’ll start the ball at a location of 20 meters from the bottom of the screen. The floor will be 4 meters from the bottom of the screen.
Let’s denote the y-coordinate of the ball after timestep t by y(t). The ball starts at timestep 0, so its initial y-coordinate is y(0). If the ball moves at constant velocity vy, then its position after timestep t is
y(t) = y(0) + t × vy ,
or in Python,
next_y = y + timestep * v_yWhat if the velocity is not constant? Then this equation won’t work. Imagine that the ball starts at velocity 0, but it is accelerating downward. At the end of the first timestep, the ball would not have moved if v_y started at 0, and that can’t be correct! However, an observation made by the mathematician Euler is that if the timestep is very small, the error won’t be very large. We’ll come back to this idea soon.
For a ball with velocity vy, with a constant acceleration ay, and timestep t,
vy(t) = vy(0) + t × ay ,
or in Python
next_v_y = v_y + timestep * a_y For our problem, ay is the acceleration due to Earth’s gravity, –9.8 meters/second2.
So here’s the idea:
In the code below, we’ve added side walls (the edges of the window, not drawn), converted all units to be relative to meters rather than pixels, and moved the velocity and position computations to their own functions.
You might wonder whether it makes much difference whether we update the position before or after we update the velocity. Here, the answer is no. Try moving the call to compute_next_velocity before the last two calls to compute_next_position. If you can notice the difference in how the program behaves, you have a keener eye than me.
The acceleration due to gravity is essentially constant near the Earth’s surface. With larger distances and long time scales, gravitational accelerations are not constant. Before looking at this problem, we’ll look at the related problem of computing the motion of a ball attached by a spring to a post.
At each timestep, we compute three quantities:
We’ll use variables x, y, v_x, v_y, a_x, and a_y to keep track of the current position and of the components of the velocity and acceleration.
Steps 1 and 3 we have already seen how to do in the previous example. Let’s look at step 2, computing the current acceleration.
The force on a body is the net acceleration on the body, multiplied by the mass of the body: f = m × a. If we knew the magnitude of the force on the ball, and the mass of the ball, we could compute the magnitude of the acceleration:
a = f/m .
For a spring, it turns out that the force is proportional to how far the string has stretched, multiplied by some constant k that depends on the particular spring. For our example, we’ll assume that the initial length of the spring was 0 (pretend it’s an infinitesmally small rubber band when not stretched). Let d be the distance of the ball from the spring attachment point, in other words, how far the spring has stretched. Now, if know that constant k, we can compute the force f:
f = k × d .
The distance d of the center of a ball at (x, y) from some point (sx, sy) is given by the good ol’ distance formula:
$d = \sqrt{(s_x - x)^2 + (s_y - y)^2}$ .
So to compute the magnitude of the acceleration a:
We’re getting there, but we have one more thing to figure out. The quantity a is just the magnitude, or size, of the acceleration. But, like velocity, acceleration has x and y components. We have to determine how much of that acceleration is in the x direction, and how much is in the y direction. The acceleration should point directly from the ball to the post where the spring is attached. But that’s not too hard. We just set the x component of the acceleration to be proportional to the ratio of the x distance of the ball from the post to the total distance:
ax = a × (sx − x)/d .
Similarly,
ay = a × (sy − y)/d .
What if we had multiple bodies, each with gravitational force exerted on it by all the other bodies? The magnitude of the force exerted on body 1 by body 2 is
$f = \displaystyle \frac{G \times m_1 \times m_2}{r^2}$,
where r is the distance between bodies in meters, m1 and m2 are the masses measured in kilograms, and G, the universal gravitational constant, is 6.67384 × 10–11.
But that’s just the force on body 1 from body 2. What about the force on body 1 from bodies 3, 4, 5, …? Because the forces are magnitudes, you can’t just add them up. But what you can do is compute the components of the acceleration of body 1 due to body 2, ax and ay, just as we did for the spring system. Then add the the components of the acceleration of body 1 due to body 3, due to body 4, body 5, etc. Once you have added accelerations due to all other bodies, you have the acceleration of body 1.
For body 2, you can compute the acceleration by computing the components of the acceleration due to bodies 1, 3, 4, 5, etc.
Once you have all the accelerations for the bodies, you can compute the next velocities for each body.