generalrelativity

This is the first in a series on simulating constraints; in this case a particle is constrained to move along a curve in planar space. To follow along, you’ll need to know a little physics, calculus and some algebra, and hopefully you’ve implemented a simple particle system before!

I’ve chosen a parabolic curve. We’re actually going to solve the problem in 2 different ways, so we’ll be seeing 2 separate representations of the curve, the first being the implicit form:

Which looks like this for x {-200…200}:

The first method we’ll use for solving the constraint is described for a circle in this paper by Andrew Witkin. There are some convenient aspects of his example, and so it may not be entirely apparent how to constrain the particle to our parabola. I’ll follow Witkin’s format to show how his formula works in this context.

First, you write a positional constraint which is just our curve equation:

The reason the curve is the constraint is because it implies that when the function’s value at (x,y) is zero, then [x,y] is located on the parabola. Next, we define the velocity constraint, which is the derivative of the positional constraint with respect to time:

If the velocity constraint is satisfied, then we’ll always end up with a valid position. Finally, we find the derivative again with respect to time. This is the acceleration constraint function:

As was the case with velocity → position, if the acceleration constraint is satisfied, so will the velocity. This means that, provided we’re starting with a valid position and velocity, we only ever need to enforce the acceleration constraint on the system.

Let’s do a quick review of Newton’s 2nd law:

Which states that force is equal to mass times acceleration. We’ll be making use of the fact that F represents the sum of all forces acting on the particle. So:

Solving for acceleration gives:

The subscripts ext and c stand for external and constraint respectively. It’s worth pointing out for clarity’s sake that when a term is bold it represents a vector, whereas italicized variables are scalars.

Now we can replace the appearances of x and y’s 2nd derivatives, breaking F into its component parts to allow the constraint force to appear explicitly:

Equation 1

The 2 unknowns now appear, but only in the 1 equation, so we can’t solve yet. If we consider a particle attached to a frictionless wire, it should be apparent that the wire, which is the source of the constraint force, can do no work and so it can not add nor remove energy from the system. The kinetic energy of the system is defined as:

To gain anything useful out of this, we have to consider what direction a force can push without changing T. From examining the derivative of this function with respect to time:

it should be apparent that the constraint force has an important relationship with the velocity vector, specifically that their dot product should equal zero (or else T would change as a direct result of the constraint force). We can use this knowledge to write the constraint force components into a 2nd equation:

Equation 2

This makes sense not just in terms of energy conservation, but also based on the constraint function chain described above. We already know our particle has a legal velocity, and so we don’t want to add any force that would ruin the circle of trust–the constraint force should only pull the external forces onto the acceleration curve.

I first solved for the x constraint component in terms of y from equation 2:

Equation 3

Solving for y’ in the velocity constraint and substituting for x’ in equation 3 reduces the variable count:

Plugging the right-hand side in for F_cx in equation 1 and solving for F_cy yields:

We could explicitly solve for F_cx too, but that’s a non-trivial number of multiplications and a divide. Luckily, we don’t have to because this is for a computational simulation; after F_cy is known, the solution to F_cx is given in equation 3.

Plugged into F=ma, we have the equations of motion for a particle constrained to the parabola x² + 100y = 0.

Parametrized Form

Now we’ll solve the problem a different way, using a reformulation of Newton’s equations called Lagrangian mechanics. Lagrangian mechanics involves finding the number of degrees of freedom (DOF) in a system, then expressing that system in terms of generalized coordinates, one for each DOF. In our example we have a single DOF, even though our particle exists in a 2-dimensional space. First, we phrase our parabola parametrically:

Equation 4

Here, q is our lone generalized coordinate. Noting that x is a vector, the transformation from the generalized coordinate to our particle’s position is given in the parametric equation (4):

The Lagrangian equations of motion for a particle experiencing only conservative forces, as is the case with a particle constrained to a wire, is:

Equation 5

David Eberly gives a full derivation of equation 5 from Newton’s 2nd law in his excellent, if often daunting, Game Physics. F, again, represents all externally applied conservative forces and T is the kinetic energy of the system, which is a function of q and q’. Before we can define T for this system, we have to derive dx/dq and x’ because kinetic energy is dependent upon velocity.

and

Now we can define the kinetic energy of the system as a function of our generalized coordinate because x’ can be expressed entirely in terms of q and its first derivative:

Having determined the kinetic energy function, we can define the rest of the terms in equation 5. We’ll see that its first term includes q’s second derivative, which is the last piece needed to advance the simulation.

Lastly, F, in this case, is just the gravitational force, which we define as the vector mg[0,1]:

Next we plug all of this back into equation 5, which then defines the Lagrangian equations of motion for a particle constrained to the parametrized parabola (q,-q²/100):

Equation 6

Solving for acceleration yields:

Now that we’ve solved for acceleration, we can advance our 2-dimensional simulation by numerically integrating our 1-dimensional differential equations. Below is a demo, showing both methods running side-by-side. If you watch long enough, you’ll see the implicit form (left) eventually leaves the parabola as error accumulates. In contrast, the parametric form will always be on the parabola, but energy will dissipate. These errors are associated with numerical integration inaccuracies and floating point imprecision.

Roll over to run simulation. Left is the implicit form, right is the parametric

Considerations

Each method has advantages and disadvantages. Both offer a relatively easy-to-follow format. Explicitly solving for a constraint force as we did with the implicit curve is more intuitive to me, but the parametric form is more computationally efficient, and the algebra isn’t as messy. A downside to the parametric form is that it isn’t a generic solution–adding new constraints to the system means solving a totally new set of equations of motion.

Below are the separate implementations that each run in a generic dynamics harness. You can see that there is more involved computationally in the implicit form:

Implicit implementation

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

var particle:Particle = Particle( p[ 0 ] );
particle.exertForce( G ); //apply gravity
 
var m:Number = particle.mass;
 
var X:Vector3D = particle.position;
var V:Vector3D = particle.velocity;
 
var Fext:Vector3D = particle.force; //generically grab all applied forces
 
var Fc:Vector3D = new Vector3D();
 
//solve for the constraint force
Fc.y = ( -50 * m * V.x * V.x - 100 * X.x * Fext.x - 2500 * Fext.y ) / ( X.x * X.x + 2500 );
Fc.x = ( Fc.y * X.x ) / 50;
 
particle.exertForce( Fc ); //exert the constraint force on the particle

Parametric implementation

1
2
3
4

var q:Number = x[ 0 ];
var qPrime:Number = v[ 0 ];
 
a[ 0 ] = ( -q * qPrime * qPrime - 50 * G * q ) / ( 2500 + q * q );

As mentioned above, both of these are prone to computational error. In the demo, I’m using an RK4 integrator. Using a Euler integration scheme results in visible error within the first second of simulation. Although RK4 does a good job of reducing error, it isn’t perfect. Ideally, an error term is introduced to keep drift in check.

drew ActionScript 3.0, Mathematics, Physics