CS184 Lecture 13 summary

Lab notes:

Its easiest to use Javascript to program the lab. In fact even if you want to use Java, you will probably find it quicker to develop via Javascript first.

Some of the machines in 330 Soda have Netscape bugs. If netscape can't load the java samples from class, you can use Explorer. To make sure Explorer finds your java classes, place them in a directory like U:\cs184-xx\myclasses and put this directory in the global classpath variable.

Dynamics

Last time we talked about kinematics as an animation tool. It allows positions (or paths) of a part of a model to be specified and the model parameters to be solved for. There is no "physics" in kinematics. Its all about shape and geometry. You can specify time trajectories with kinematics, but they're really like geometric objects (curves) in space-time.

Dynamics covers the full gamut of physical motion, including Newton's laws, elasticity, friction, motors, muscles etc. Its a big subject, and we give only the flavor of it here.

Physics 0.1

Let's recall a few facts from Newtonian mechanics. First of all:

F = m a

where F is the force on an object, m is its mass, and a is its acceleration. This law is particularly simple. Mass is a scalar, so force and acceleration are parallel. It works for rigid bodies as well as point masses.

To simulate a system with dynamics, assuming you know the applied force, you can integrate a to get the velocity v, and then integrate v to get the position x. In a real simulator, you would use discrete time steps and approximate integrals with sums like v' = v + a dt, and x' = x + v dt.

Numerical integration has problems however, in particular error accumulation. Adding up many small numbers means adding up many small error terms, which can become a big error term. Real physical systems satisfy conservation laws. For translation, we have:

Conservation of Energy: in the absence of force, kinetic energy K is conserved:

K = 1/2 m v2

Conservation of Momentum: in the absence of force, momentum p is conserved:

p = m v

So a good starting point for a numerical scheme is that it satisfy conservation laws. It doesnt matter so much for blocks on a table, but it does for spacecraft or satellite simulation.

Closed-form Dynamics

A handful of dynamical systems have an easily-computable closed-form solution. (most dont have a closed-form solution at all). Simple spring-mass and pendulum systems can be solved. Recall that a spring applies a force proportional to its displacement, so

F = - k x

Since F = ma and a = d2x/dt2, we have an ordinary differential equation:

d2x/dt2 = - k/m x

which has a solution with sinusoids x = A sin(w t + j)  where w2 = k/m.

Here is an example environment with a mass-spring system. Its actually kinematically-generated motion with a script.

Rigid-Body Dynamics: Rotation

The story for rotational dynamics is more complicated than for translation, and there are very complex and interesting behaviors. Think about spinning tops, footballs flying through the air, riding a bicycle (with or without holding the handlebars). Instead of Newton's laws, for rotation we have Euler's equation:

Ia + w x I w = T

where T is an external applied torque (or moment), w is the angular velocity, a is its derivative (the angular acceleration) and I is the inertia matrix, which plays the role of mass. Unlike mass, which is a scalar, the inertia matrix is a 3x3 symmetric matrix. It depends on the shape of an object, and has different values in different coordinate systems.

As for translation, we have conservation laws for rotation:

Conservation of angular momentum: In the absence of a  torque T, angular momentum p is conserved:

p = I w

Conservation of energy: In the absence of a torque T, kinetic energy K is conserved:

K = 1/2 wT I w

To simulate a rigid body, we could start with Euler's equation. Assuming T is known, at each time step, we would compute the angular acceleration a. Since a is the time derivative of w, we can integrate it to get the change in w. But how do we "integrate" w?

Recall that angular velocity w represents rotation about an axis parallel to w, and with angular speed (radians per second) proportional to the magnitude of w. We want to make a small change in an object's orientation as given by an R matrix. Or equivalently, we want to know dR/dt as a function of w.

To do this, we use the axis-angle formula to get an approximation for a small rotation by w dt. The rotation matrix is approximately.

I + (wx) dt

To get dR/dt we apply this small rotation to the current orientation R, then subtract of the old value of R, and divide by dt.

dR/dt = ((I + (wx) dt)R - R) / dt

so the time derivative of R is

dR/dt = (wx)R

To compute the changes in R with time, we could integrate dR/dt. Because of numerical errors, it works better to turn (w dt) into a real rotation matrix and multiply it by R. Some sample code that does this follows soon.

Symplectic Integrators

Integrating the Euler equations will probably violate both conservation of momentum and energy after some amount of time. Its also a complicated equation. There is a simpler and more accurate approach based on symplectic integration. You dont need to worry about what it means, just that it guarantees conservation by using the conservation equations as constraints. We have conservation of momentum:

I w = p

But both w and I are varying with time. However, in the object's coordinate frame, the inertia matrix I0 is constant. In this coordinate frame, the velocity is w' = RT w, and p' = RT p. So we have

I0 RT w =  RT p

and rearranging gives

w = R I0-1 RT p

All these quantities are known, so we can compute w. Then we can use the expression above for dR/dt to update R to the next position. This simple scheme is guaranteed to conserve momentum.

Unfortunately, it is still possible to accumulate errors that lead to large energy changes. Here is a dynamic simulation using the above scheme. The right-hand bar shows the (normalized) energy of the system over time. Note the dramatic increase in energy as the block spirals in.

The energy of the system is 

1/2 wT I w  =  1/2 pT w

by making a small change to the orientation R at each step to correct fluctuations in this energy term, we can guarantee that energy is conserved. The correction is an infinitesimal "spin" normal to both p and w. Here is the same environment, with an energy correction term.

The Java source code for the dynamics integrator is here. It uses some classes for rotations and basic matrix operations in java that you may find useful. Click on these links to grab that code.

Chaotic behavior

Many simple systems, like multi-link pendulums, dont have a closed-form solution. Many dynamical systems including those exhibit chaotic behavior. Chaotic systems have the property that their trajectories are aperiodic (not repeating) and their state is highly sensitive to the initial conditions. So a small change to the initial conditions of a chaotic system lead to large changes in the outcome over time.

For such systems, its  a good idea to run simulations with small random changes to the model. Because simulations are deterministic, you get the same outcome every time you run the simulation. If you ran a real chaotic system several times, you would almost certainly get very different outcomes each time. Video examples.