Here is the robotics approach to this problem.
What you need is the kinematic relationship between bodies that always satisfies the joint. This relationship connects the joint coordinates (angles, distances, etc) and their derivatives to the cartesian motion vectors of each body. The goal is to define the configuration of the system in terms of the joint degrees of freedom.
Then you can reduce the system into an ODE in terms of the joint coordinates which you integrate. As a post processing step for each time frame $t$, joint positions $\boldsymbol{q}$, joint velocities $\boldsymbol{\dot q}$ and joint accelerations $\boldsymbol{\ddot q}$ can yield the full motion and force description of each body.
The easiest problems to formulate are those of where each body is connected with a joint to a single parent body, and the common parent of all bodies is the immovable ground. These tree like systems have recursive algorithms that produce almost analytical results.
Next in complexity would be problems where there are kinematic loops joining different branches of the tree together that require additional constraint forces to enforce all the loops.
Finally is the problem of a free floating system of bodies which may not have a common parent body to define all other body configurations relative to. In this case you have to mix the robotic joint coordinates $\boldsymbol{q}$ with a standard rigid body parametrization of position and Euler angles, or quaternions and integrate everything at the same time.