Date: Mon, 16 Dec 1996 23:40:02 GMT Server: NCSA/1.5 Content-type: text/html Last-modified: Wed, 28 Feb 1996 13:45:23 GMT Content-length: 14462 dynamics

Using forces to animate objects


Introduction

There are several good reasons why you may want to compute motion of objects based on forces, rather than attempt a kinematic description of object motion.


Background

The basis of many particle systems is a description of the acceleration of the particles. If the accelerations are given, then the system is called second order because two integrations must be done to get the positions.

As an example, consider a mass on the end of a spring sliding on a frictionless horizontal table. the force on the mass is given by:

F = -k*x

where F is the force, k is the spring constant and x is the length of the spring. The negative sign appears because the force from a spring pulls in the opposite direction to the spring extension.

Since F = m*a where a is the acceleration,

a = -(k/m)*x

Now given the acceleration, and if we know an initial velocity and position, we can find the position of the mass at any later time. In the case of this linear spring, we can integrate the acceleration directly to find that the position is a sinusoidal function of time. The constant (k/m) will determine the frequency of the sinusoidal function, while the initial position, x0, and velocity, v0, will determine the amplitude and phase [5]. The position is given by:

In this simple case, the solution is available, so one could code a kinematic solution (that is, just plug time into the sin function), but notice that the description of the system in differential form is more compact and straightforward. Also if we changed the force to a more realistic "stiff spring" model such as:

F = -k*x - c*x^3

we could not directly integrate it, but would be forced to use numerical methods. The next section will consider how to perform the integration.


Numerical integration

Solving a dynamic system means figuring out how to move the system forward in time from some set of initial conditions to compute the positions as a function of time. For instance you might want of trace the trajectory of a ball after it is fired from a cannon. In general, direct analytical solutions of the differential equations governing a system are hard or impossible. We will outline some techniques for solving a system by simple (but perhaps not optimal) numerical methods.

As a disclaimer, it is difficult to construct general numerical schemes that work for various physical systems. What follows is a specific approach which works for some systems. Consider it the barest introduction to the topic of numerical introduction. Much more detail may be obtained in [1] or [3] or many other books.

We will consider three specific second order systems.

For each system we will first introduce the physics (the force law) then the scheme to integrate the force law to produce the time-dependent positions. We will start with a integration scheme which will be used only for intrducing the notion of numerical integration. While useful to look at, it should almost never be used for actual work.
The Euler integration algorithm

Let me state again that this algorithm, while easy to understand, is so inefficient that it should only be used whem programming effort dominates throughput. We want to start with a differential equation and show how to integrate in numerically. Start with the differential equation

We want to solve this equation by stepping time discretely forward in small steps ("small" still needs to be defined). One discrete approximation of this equation is:

where the n and n-1 refer to two time steps sufficiently close together that the differences approximate the derivitive. Rearranging this equation yields:

which gives a explicit form for stepping the system from time n-1 to time n. This form is thus a way to step a dynamic system forward in time given an initial state x0. The inefficiency of this method occurs because of the assumption of constant f(x) across the entire time step. To be accurate the time step has to be very short compared to the natural time constants of the system.


The Verlet integration algorithm

A much better alogrithm can be derived by averaging slopes across the small time interval of integration. One such formulation is called the velocity form of the Verlet algorithm [1]. This algorithm is most useful for second order systems where the force on an object is specified and the position is desired as a function of time. Since the accelerations are specified, one integration must be done to calculate the velocity and a second done to calculate the position.

The Verlet scheme is first updates the position, then use the old and new position information to update the velocity.

Where are the accelerations calculated from a force law describing the system and which are usually a function of position. Note that an is calculated from the just-completed calculation for xn. Evaluating both of these equations advances the system one time step. To start the solution, a velocity and position have to be known at t=0. Choosing delta t requires some care. Generally, you can start with a time step which is around .2 of the fastest time-constant in the system.

The two Verlet equations can be vector equations if the motion is 2 or 3 dimensional.

Now we are ready to consider some specific physical systems.


Gravitation and the three-body problem

The accleration due to a gravitating mass is

where G is a strength contant, M is the mass of the object pulling on you, and r is the distance between you and M. Since a and r are actually vectors we must derive a form of this equation which is useful in 2 or 3 dimensions. We will use 2D here so we need equations for the x and y components of the acceleration. Assuming body one is located at r1 and body two is located at r2, and that we want the acceleration of body one:

Where theta is the angle measured from the x-axis of the vector r pointing toward body one from body two.

Converting to a Cartesian form the acceleration of body one is:

where the vectors

The calculation procedure for each time step is to compute:

  1. the acceleration, a, based on the positions at time n-1
  2. a new set of positions using the Verlet method
  3. the acceleration at time n based on the newly computed positions
  4. a new velocity from the Verlet method using the acclerations at times n-1 and n.
For the gravitational animation given at the beginning of this page, three masses were simulated. The accleration of each mass was determined as the vector sum of the accelerations caused by each of the other two bodies.
Water Waves

The water wave solver presented here is based on the derivation in Kass and Miller [2] for the case of shallow water, low amplitude, non-breaking waves. The form of the resulting equation of motion for the waves looks like the classical linear wave equation with a propagation velocity proportional to the depth of the water.

where h is the height of the water surface and d is the depth, that is, d(x,y) = h(x,y)-b(x,y) where b is the vertical position of the containing vessel (or ocean bottom) at (x,y). The costant g is poroportional to the force of gravity. This partial differential equation can be discretized in many ways. We will use a method which is stable enough to be acceptable for computer graphics.

First note that the partial with respect to time is the acceleration of a small surface element of water, so that if the right side of the equation can be put in a discrete form we can apply the Verlet method to a 2D grid of "bodies", each representing a small chunk of water.

To solve this equation numerically we need to discretize it both in space, on a two dimensional grid, and in time. The discrete spatial approximation is that the second derivitives (at time n) are

were i and j are the grid indices in the x and y directions and n is the time index. Note that at the edges of the array the discrete partial derivitives depend on values outside the array. These boundary conditions need to be specified for a solution. One easy boundary condition, corresponding to no transport across the coundary, is to copy the value at the edge of the array whenever a value outside the array is needed.

It seems that the solution is more stable if the minimum depth of the water is limited to about .001 of the average value, rather than letting it go to zero. So the depth at time n is given by

To get the surface motion, first assign an initial height and vertical velocity to each grid point (i,j) then at each time step compute :

  1. the acceleration, a, based on the water heights at time n-1

  2. a new set of positions using the Verlet method for each (i,j)
  3. the acceleration at time n based on the newly computed heights

  4. a new velocity from the Verlet method for each (i,j)
  5. Ensure conservation of volume of the total volume by adjusting the average water height to a constant. That is, sum all the heights over the whole grid, then correct the sum to equal the initial sum by adding a small increment to each grid location (i,j).
These steps are repeated to generate moving waves.

Disclaimer: While this integration method seems to give visually reasonable results, it should not be used for analytical simulations without careful validation. The wave equation is notoriously hard to integrate using explicit techniques, like the one described here.


Billards

A billards, hard ball, system is different than the other two systems just described because the forces between balls are zero until they touch and become very large if the balls try to pass through each other. This means that there is a large, impulsive force just when the balls meet and at no other times. The Verlet integration scheme will fail badly on this system because the accelerations are large for a short time and thus not smooth enough to average. What we will do to step the billards system forward in time is to calculate the total change in velocity from a collision, without worrying exactly how forces change the velocity.

The method described here is a less exact version of that described in [4]. The method described in the paper steps time at uneven intervals, and is thus unsuitable for animation. My modification steps time uniformly, at the expense of exact collision dynamics.

The change in velocity during impact can be derived for frictionless balls of equal mass by noting that the the impact force must act in a direction parallel to the line connecting the centers of the two impacting balls. The change in velocity must be parallel to the connecting line also, with the velocity component parallel to the line having its sign reversed by the collision and the velocity component perpendicular to the line unchanged. Projecting the initial velocity onto the line connecting the centers, negating the result, and resolving it back into x and y velocity components gives the velocity change. If i and j are the indices of the colliding balls define:

then delta v for ball i is given by the following where the right-most term represents the projection of the velocity onto the line and the other term converts the projection back to x,y coordinates.

The calculation procedure for each time step is to compute:

  1. the delta v based on the positions of the balls
  2. a new set of positions using

  3. if walls are present detect collisions and modify velocities.
The time step needs to be small enough so that the balls do not penetrate each other too much during one time step.
References

  1. An Introduction to Computer Simulation Methods, Part 1, by Harvey Gould and Jan Tobochnik, Addison-Wesley, 1988

  2. Rapid, Stable Fluid dynamics for Computer Graphics, by Micheal Kass and Gavin Miller, Computer Graphics, Vol 24 #4, Aug 1990, pp 49-55

  3. Numerical Recipes, by William Press, Saul Teukolsky, William Vetterling and Brian Flannery, Cambrige, 2nd edition, 1992

  4. Studies in Molecular Dynamics. I. General Method, by B. Alder and T. Wainwright, Journal of Chemical Physics, Vol 31 #2, Aug 1959, pp 459-466

  5. Physics for students of science and engineering, Part 1, by Robert Resnick and David Halliday, Wiley, 1963.


Comments about Theory Center online documents are welcome and may be sent to doc-comments@tc.cornell.edu.

Last modified, 10/12/95 B. Land.

Copyright Statement