A novel symbolic formulation is presented to model dynamics of large-scale open-loop holonomic multibody systems, by using absolute joint coordinates and via matrix transformation, instead of solving constraint equations. The resulting minimal set of second-order linear ordinary differential equations (ODEs) can be used for linear vibration analysis and control directly. The ODEs are generated in three steps. Firstly, a set of linearized ODEs are formulated in terms of absolute coordinates without considering any constraint. Secondly, an overall transform matrix representing constraint topology for the entire constrained system is generated. Finally, matrices for a minimal set of ODEs for the open-loop holonomic multibody system are obtained via matrix transformation. The correctness and efficiency of the presented algorithm are verified by numerical experiments on various cases of holonomic multibody systems with different open-loop topologies, including chain topology and tree topology. It is indicated that the proposed method can significantly improve efficiency without losing computational accuracy.