In the limit of infinite stiffness, the differential equations of motion of stiff mechanical systems become differential algebraic equations whose solutions stay in a constraint submanifold $\widehat{\mathcal{P}}$ of the phase space. Even though solutions of the stiff differential equations are typically oscillatory with large frequency, there exists a slow manifold $\widetilde{\mathcal{P}}$ consisting of nonoscillatory solutions; $\widetilde{\mathcal{P}}$ has the same dimension as $\widehat{\mathcal{P}}$ and converges to it as the stiffness approaches infinity. We introduce an iterative projection algorithm, IPA, that projects points in the phase space of a stiff mechanical system onto the associated slow manifold $\widetilde{\mathcal{P}}$. The algorithm is based on ideas such as micro-integration and filtering coming from the field of multiscale simulation and is applicable to initializing integration algorithms for both stiff ODEs and DAEs, including the initialization of Lagrange multipliers. We also ...