Abstract: A mixed variable formulation and the relevant precise integration method are proposed for the dynamic response analysis of multilayered road structures. By performing the Fourier-Bessel transform, the partial differential wave motion equation in the frequency wave number domain can be decoupled into two sets of second-order ordinary differential equations, one for P-SV components and the other for SH components. By introducing the dual vectors of stress and displacement, the second-order ordinary differential equation is further reduced to a homogeneous first-order one. The solution is in the form of an exponential function. By employing the precise integration method, very high accuracy can be achieved. Furthermore, the mixed variable formulation of the solution of wave motion equation facilitates the assembly of layers and improves the computational efficiency. The proposed method is applicable to arbitrary distribution of loads. The computation is stable and convenient for the computer programming. The accuracy and rationality of the proposed method are verified by comparing the solutions with the BISAR software and the experimental results. Some numerical results are presented to reflect the deformation characteristics and stress distribution of the road structures under the tire loads. The research results are useful for the design of road pavement structures.