This research proposes an implementation method of linearized equations of motion for multibody systems with closed loops. The null space of the constraint Jacobian is first pre-multiplied to the equations of motion to eliminate the Lagrange multiplier and the equations of motion are reduced down to a minimum set of ordinary differential equations. The resulting differential equations are functions of ail relative coordinates, velocities, and accelerations. Since the coordinates, velocities, and accelerations are tightly coupled by the position, velocity, and acceleration level constraints, direct substitution of the relationships among these variables yields very complicated equations to be implemented. As a consequence, the reduced equations of motion are perturbed with respect to the variations of all coordinates, velocities, and accelerations, which are coupled by the constraints. The position, velocity and acceleration level constraints are also perturbed to obtain the relationships between the variations of all relative coordinates, velocities, and accelerations and variations of the independent ones. The perturbed constraint equations are then simultaneously solved for variations of all coordinates, velocities, and accelerations only in terms of the variations of the independent coordinates, velocities, and accelerations. Finally, the relationships between the variations of all coordinates, velocities, accelerations and these of the independent ones are substituted into the variational equations of motion to obtain the linearized equations of motion only in terms of the independent coordinate, velocity, and acceleration variations.