Abstract: According to the precise integral method, a new numerical algorithm is proposed for dynamical responses of viscoelasticity structure with the fractional derivative constitutive relation. In this approach, the fractional differential equation for the dynamic of a system is transformed into a set of first order ordinary differential equations which contain fractional derivative terms. The precise integral method is used to integrate these terms and obtain the response of the system. Numerical results obtained using this scheme agree well with those obtained using an analytical method and an numerical scheme proposed by ZhangShimizu. Results also show that the solution converges as the step length is decreased.