In this paper, we present efficient algorithms for computation of the residual of the constrained discrete Euler–Lagrange (DEL) equations of motion for tree structured, rigid multibody systems. In particular, we present new recursive formulas for computing partial derivatives of the kinetic energy. This enables us to solve the inverse dynamics problem of the discrete system with linear computational complexity. The resulting algorithms are easy to implement and can naturally be applied to a very broad class of multibody systems by imposing constraints on the coordinates by means of Lagrange multipliers. A comparison is made with an existing software package, which shows a drastic improvement in computational efficiency. Our interest in inverse dynamics is primarily to apply direct transcription optimal control methods to multibody systems. As an example application, we present a digital human motion planning problem, which we solve using the proposed method. Furthermore, we present detailed descriptions of several common joints, in particular singularity-free models of the spherical joint and the rigid body joint, using the Lie groups of unit quaternions and unit dual quaternions, respectively.