Traditionally, flexible multibody dynamics problems are formulated as initial value problems: initial states of the system are given and solving for the equations of motion yields the dynamic response. Many practical problems, however, are boundary rather than initial value problems; two-point and periodic boundary problems, in particular, are quite common. For instance, the trajectory optimization of robotic arms and spacecrafts is formulated as a two-point boundary value problem; determination of the periodic dynamic response of helicopter and wind turbine blades is formulated as a periodic boundary value problem; the analysis of the stability of these periodic solutions is another important of problem. The objective of this thesis is to develop a unified solution procedure for both initial and boundary value problems. Galerkin methods provide a suitable framework for the development of such solvers. Galerkin methods require interpolation schemes that approximate the unknown rigid-body motion fields. Novel interpolation schemes for rigid-body motions are proposed based on minimization of eighted distance measures of rigid-body motions. Based on the proposed interpolation schemes, a unified continuous/discontinuous Galerkin solver is developed for the formulation of geometrically exact beams, for the determination of solutions of initial and periodic boundary value problems, for the stability analysis of periodic solutions, and for the optimal control/optimization problems of flexible multibody systems.