In-class exercise: simulating gravity between two point masses

You've succeeded in following the motion of a skydiver as he falls towards Earth. That was pretty simple: the gravitational force on the diver was always exactly the same:



      grav force  F  =  m * g  =  m * (-9.8 m/s^2)

Sure, if you want to get technical, the gravitational force decreases slightly with height; but even at an altitude of 10,000 meters, it doesn't make much of a difference.

But consider a different situation: two bodies, alone in the universe. They interact only through their gravitational force on each other. In this case, the size of gravitational force is


                         G * m1 * m2 
      grav force  F  =  -------------
                            r * r

where


  G  =  universal constant of gravitation  =  6.67 x 10^(-11) N*(kg*kg)/(m*m)

  m1 =  mass of first object (kg)

  m2 =  mass of second object (kg)

   r =  distance between objects (m)

This may seem like an artificial situation -- two objects alone in the universe? -- but it's actually quite a decent approximation to the real Solar System: the gravitational force between the Sun and the Earth is much, much larger than that between the Earth and any other planet. The dynamics of celestial systems are ususally dominated by gravity, and, because the distances between objects are usually much, much larger than the sizes of the objects themselves, we can treat them all as point masses. Or can we?

Let's find out. Suppose that there are two objects, initially at rest, a distance of 100 meters apart from each other.

To simplify things, we fix the massive object (number 1) in place at x = 100 m, and allow only the tiny object (number 2) to move. Please write down on paper your answers to the following questions:



  Q: What will happen when we release the tiny object?

  Q: Suppose that both objects are really "point masses", so that 
         they can pass each other without touching.  How far to the
         right will the tiny object go before it ought to stop
         and come back again?

  Q: What is the initial Kinetic Energy of the tiny object?

  Q: What is the initial Gravitational Potential Energy of the tiny object?

Now, let's simulate the motion of the tiny object. Take a look at this half-written program:

inclass_grav.m


function inclass_grav(timestep, duration)
%
% A very simple illustration of the dangers of numerical integration
%   of gravity.  We follow the motion of a very massive object
%   (fixed in place) and a tiny object (free to move).  They start
%   at rest.  The tiny object falls towards the big one.
%   We follow the motion for "duration" seconds with the given "timestep"
%
%           Arguments:    timestep   (input)       size of timestep (seconds)
%
%                         duration   (input)       number of seconds to 
%                                                     run the simulation
%
% MWR 4/28/2003

% mass of each object (kg)
m1 = 1e9;
m2 = 1;

% initial position of each object (m).  Note that
%   the massive object is fixed in place and never moves
x1 = 100;
x2 = 0;

% initial velocity of each object (m/s).  Again, note that the
%   massive object is fixed in place and can never move
v1 = 0.0;
v2 = 0.0;

% gravitational constant, in MKS units
G = 6.67E-11;

% initial total energy: GPE + KE
dist = x2 - x1;
initial_energy = -G*m1*m2/abs(dist) + 0.5*m2*v2*v2;
fprintf(1, 'initial energy is %9.4e\n', initial_energy);

t = 0;
while (t < duration) 

  % calculate the distance between the object here ...
  dist = 

  % calculate current total energy
  current_energy = 

  % now find the gravitational force between the objects ...
  force = 

  % and the acceleration of the tiny object
  accel2 = 

  % use a simple Euler method to calculate new velocity and position
  new_v2 = 
  new_x2 = 

  % print out information on the CURRENT (not new) quantities
  fprintf(1, ' %9.4f   %10.5f  %10.5f  %9.4e \n', t, x2, v2, current_energy);

  % update variables for next time through the loop
  v2 = new_v2;
  x2 = new_x2;
  t = t + timestep;
end

% end of program

It employs Euler's method to calculate the position and velocity of the tiny object as a function of time. But it's missing some pieces in its inner section. Your job in class today:

Copy this program to your desktop
Add the pieces which are missing
Run the program with the following input arguments:
- timestep = 100, duration = 1000
- timestep = 100, duration = 2000
- timestep = 100, duration = 3000
- timestep = 100, duration = 4000
- timestep = 100, duration = 5000
- timestep = 100, duration = 10000
Each time you run the program, compare the position, velocity and total energy of the tiny object to your predictions above.