We can model the movement of an isolated object using $F = m\ \ddot{x}$ and $\tau = I \dot\omega + \omega \times (I\ \omega)$.
But objects are often connected and are constrained to move together.
Here we'll look at a particular constraint and how to simulate it.
Derivatives of local vectors
First we'll look at the derivatives of vectors attached to moving objects, as these are necessary in formulating constraints.
Below is an object rotating with angular velocity $\omega$ about its centre of mass, $x$. A point $c$, on the surface of the object, is defined by the vector from the origin of the local coordinate system (i.e. the centre of mass) to the point on the boundary.
The point, $c$, has velocity $\dot{c}$: $$\begin{array}{rcl} \dot{c} &=& {d \over dt}\ c \\ &=& w \times c \end{array}$$
Its acceleration is $\ddot{c}$: $$\begin{array}{rcl} \ddot{c} &=& {d \over dt}\ \dot{c} \\ &=& {d \over dt}\ (\omega \times c) \\ &=& \left({d \over dt}\ \omega\right) \times c + w \times \left( {d \over dt}\ c \right) \\ &=& \underbrace{\dot\omega \times c}_\textrm{tangential} + \underbrace{\omega \times (\omega \times c)}_\textrm{centripetal} \end{array}$$
Furthermore, $${d \over dt}\ (x+c) = \dot{x} + \omega \times c$$
and $${d^2 \over dt^2}\ (x+c) = \ddot{x} + \dot\omega \times c + \omega \times (\omega \times c)$$
Point-to-point constraint
A point-to-point constraint ensures that a point, $c_1$, on one object remains in contact with a point, $c_2$, on another object. This can be considered a "ball in socket" or "hip joint".
To maintain contact between $c_1$ and $c_2$, we need a force acting at that point of contact. The force acts positively on one object (say $O1$) and negatively on the other object (say $O2$):
We must include this force in the equations of motion of the two objects. Note that the force induces a torque, $c \times F$. We'll also include a gravitational force, $m\ g$: $$\begin{array}{rcl} F + m_1\ g &=& m_1\ \ddot{x_1} \\ c_1 \times F &=& I_1\ \dot\omega_1 + \omega_1 \times (I_1\ \omega_1) \\ \\ -F + m_2\ g &=& m_2\ \ddot{x_2} \\ c_2 \times (-F) &=& I_2\ \dot\omega_2 + \omega_2 \times (I_2\ \omega_2) \\ \end{array}$$
Each equation above is a vectorial equation, so is really three scalar equations in $x$, $y$, and $z$.
Now we have 15 scalar unknowns $$\ddot{x_1}, \dot\omega_1, \ddot{x_2}, \dot\omega_2, F$$
and only 12 scalar equations (3 for each equation above).
We need three more equations if we are to have a complete system of equations, which we can use to solve for the unknowns.
So we will add a point-to-point constraint, which requires that the accelerations of the $c_1$ and $c_2$ points be identical. $${d^2 \over dt^2}\ (x_1+c_1) = {d^2 \over dt^2}\ (x_2+c_2)$$
But the constraint has to be in the unknowns of the system of equations, so we will instead say $$\ddot{x_1} + \dot\omega_1 \times c_1 + \omega_1 \times (\omega_1 \times c_1) = \ddot{x_2} + \dot\omega_2 \times c_2 + \omega_2 \times (\omega_2 \times c_2)$$
Now the system of equations is $$ \left[\begin{array}{ccccc} m_1 & 0 & 0 & 0 & -1 \\ 0 & I_1 & 0 & 0 & -\widetilde{c_1} \\ 0 & 0 & m_2 & 0 & 1 \\ 0 & 0 & 0 & I_2 & \widetilde{c_2} \\ 1 & -\widetilde{c_1} & -1 & \widetilde{c_2} & 0 \end{array}\right] \ \left[\begin{array}{c} \ddot{x_1} \\ \dot\omega_1 \\ \ddot{x_2} \\ \dot\omega_2 \\ F \end{array}\right] = \left[\begin{array}{c} m_1\ g \\ - \omega_1 \times (I_1\ \omega_1) \\ m_2\ g \\ - \omega_2 \times (I_2\ \omega_2) \\ - \omega_1 \times (\omega_1 \times c_1) + \omega_2 \times (\omega_2 \times c_2) \end{array}\right] $$
where each entry in the left-hand matrix is, itself, a $3 \times 3$ matrix. For example, $$m_i = \left[\begin{array}{ccc} m_i & 0 & 0 \\ 0 & m_i & 0 \\ 0 & 0 & m_i \end{array}\right]$$
and $\widetilde{c_\textrm{ }}$ is a skew-symmetric matrix $$\widetilde{c_\textrm{ }} = \left[\begin{array}{ccc} 0 & -c_z & c_y \\ c_z & 0 & -c_x \\ -c_y & c_x & 0 \end{array}\right]$$
which has the property that $\widetilde{c_\textrm{ }}\ \omega = c \times \omega$.
We can now solve this $15 \times 15$ system to get values for $$\ddot{x_1}, \dot\omega_1, \ddot{x_2}, \dot\omega_2, F$$
These are the accelerations that go in to the state derivative, $\dot{y}$, from which we can compute $y(t + \Delta t)$.
Then the objects will move while satisfying the constraint.
Other constraints
Other constraints are done similarly: Add forces and torques to physically maintain the constraint, and add constraint equations to define the constraint.