The premise of this work is that the presence of high stiffness and/or frictional contact/impact phenomena limits the effective use of high order integration formulas when numerically investigating the time evolution of real-life mechanical systems. Producing a numerical solution relies most often on low-order integration formulas of which the paper investigates three alternatives: Newmark, HHT, and order 2 BDFs. Using these methods, a first set of three algorithms is obtained as the outcome of a direct index-3 discretization approach that considers the equations of motion of a multibody system along with the position kinematic constraints. The second batch of three algorithms draws on the HHT and BDF integration formulas and considers, in addition to the equations of motion, both the position and velocity kinematic constraint equations. Numerical experiments are carried out to compare the algorithms in terms of several metrics: (a) order of convergence, (b) energy preservation, (c) velocity kinematic constraint drift, and (d) efficiency. The numerical experiments draw on a set of three mechanical systems: a rigid slider-crank, a slider-crank with a flexible body, and a seven body mechanism. The algorithms investigated show good performance in relation to the asymptotic behavior of the integration error and, with one exception, result in comparable CPU simulation times with a small premium being paid for enforcing the velocity kinematic constraints.