In this thesis, we propose several algorithms for model reduction of second order dynamical systems. These various projection methods are based on singular value decomposition, Krylov projection, and balanced truncation. In many cases models are given in second order form, and the goal is to produce a reduced order system which is in second order form, gives an accurate approximation of the original system, and maintains some important properties such as stability and passivity. Model reduction on first order linear time invariant dynamical systems has been extensively studied, algorithms and theory are well-developed. People usually study and deal with second order system by transforming to first order form which doubles the dimension. This can be inefficient and generally does not respect the second order form. The reduced model is not realizable as a second order system. So far only a very few algorithms have been proposed for second order model reduction. Most of these are not practical for large scale settings, and no error bounds have been provided. In this thesis, a global error bound is given for some of the algorithms based on SVD and balanced truncation, the error bound is bounded by a constant times the summation of the neglected singular (or Hankel singular) values, that means those second order model reduction algorithms provide accurate approximations to the original systems. The structures of controllability and observability Gramians P and Q are discussed. All algorithms developed in this thesis have been implemented and shown to be numerically efficient, and applicable to large scale settings. All algorithms are implemented in Matlab, some of them are implemented in Fortran and C separately for which we use LAPACK. In this thesis, we apply our algorithms to three real models. The performance of our algorithms is compared with some of the previously existing algorithms. It turns out that most of our algorithms are very competitive with existing methods.