This investigation is the first contribution of a two-part research work concerning the theoretical development of a multibody approach to analyze the constrained dynamics of articulated mechanical systems. In this paper, a method for investigating the linear and nonlinear stability of the dynamic behavior of mechanical systems modeled as multibody systems subjected to holonomic and nonholonomic constraints is presented. To this end, the nonlinear equations of motions that assume a complex index-three differential-algebraic form are systematically formulated and directly linearized by using an automatic procedure based on a hybrid symbolic-numeric approach devised in this work. The proposed stability analysis method, therefore, is based on the formulation of a generalized eigenvalue problem and represents a viable computer-aided approach suitable for analyzing multibody mechanical systems having different degrees of complexity. Furthermore, an extension of the generalized coordinate partitioning algorithm is introduced in this paper for handling nonholonomic multibody systems leading to a robust and general multibody computational procedure referred to as the Robust Generalized Coordinate Partitioning Algorithm (RGCPA). Since the methodologies employed in this paper to study the stability of multibody mechanical systems are general and versatile, they can be easily implemented in general-purpose multibody computer programs and readily used to analyze several mechanical applications having engineering interest.