The discrete variational method based on Hermite interpolation is studied for the simulation problem of multibody system dynamics. Hermite interpolation is carried out to interpolate the generalized coordinate and the generalized velocity, firstly. Then, using the Gauss numerical integral method, the Hamilton principle and discrete variational principles of mechanics, the discrete variational mathematical model based on the Hermite interpolation containing known and unknown derivative information is established. High accuracy simulation results is obtained by solving the presented discrete variational equations. This method can accurately maintain the constraint equations on the larger time step, and keep the total energy of the system in a certain range of bounded variation. Therefore, it can be applied for the long time simulation.