In the numerical integration of ordinary differentialequations (ODEs) in structural dynamics community, \theta _1 method hascharacteristics of controlled numerical dissipation and second-orderaccuracy for systems with or without physical damping. Based on thesecharacteristics, \theta _1 method is extended to the numerical integrationof motion equations in multibody system dynamics. The solved motionequations are index-3 differential-algebraic equations (DAEs) andindex-2 over-determined DAEs (ODAEs). Numerical experiments validate the\theta_1 method, experiments also show the relationship of numericaldissipation with parameter \theta_1.As for theintegration of index-3 DAEs by \theta _1 method, it has violation ofvelocity constraint, while for index-2 ODAEs, there are no violation ofposition and velocity constraint in the view of computer precision. Inaddition, experiments illustrate that, for non-conservative system motionequations in the form of index-3 DAEs and index-2 ODAEs, \theta _1 methodhas second-order accuracy. In the end, \theta _1 methods for motionequations are compared with other direct-time integrations from the CPU timepoint of view.