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 all relative coordinates, velocities, and accelerations. Since the variables are tightly coupled by the position, velocity, and acceleration level coordinates, 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 variables, 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 variables only in terms of the variations of the independent variables. Finally, the relationships between the variations of all variables 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 variables variations.