This paper presents an analytical/numerical method for linearizing the equations of motion and evaluating the system Jacobian matrices of mechanical systems with closed chains. The linearization algorithm developed here first identifies and linearizes basic recursive kinematic relationships and then applies the chain rule to the derivation of the equations of motion under the framework of recursive formulation. This method can be incorporated into formulating recursive equations of motion for general multibody dynamic systems, to handle large scale systems. Since no numerical differentiation is used in the proposed algorithm, its accuracy is comparable to symbolic, closed-form linearization. Moreover, without the need of repetitious computation to select proper perturbation quantities, this method is computationally more efficient than the finite difference method.