The main aim of this paper is to compare the effectiveness of numerical techniques used for spatial multibody dynamics simulation by applying the natural complement method. In the paper, seven numerical schemes are considered: zero eigenvalue formulation, Pseudo Upper Triangular Decomposition, Schur decomposition, Singular Value Decomposition, QR decomposition, coordinate partitioning and Wang-Huston formulation. In order to illustrate the effectiveness of these schemes, two McPherson struts are considered. Simulations are performed with four error tolerance values and for three stabilization cases. Some suggestions on possible applications of the selected methods are formulated.