The description of real-life engineering structural systems is usually associated with some amount of uncertainty in specifying material properties, geometric parameters and boundary conditions. In the context of structural dynamics it is necessary to consider joint probability distribution of the natural frequencies in order to account for these uncertainties. Current methods to deal with such problems are dominated by approximate perturbation methods. In this paper a new approach based on an asymptotic approximation of multidimensional integrals is proposed. A closed-form expression for general order joint moments of arbitrary number of natural frequencies of linear stochastic systems is derived. The proposed method does not employ the small-randomness and Gaussian random variable assumption usually used in the perturbation based methods. Joint distributions of the natural frequencies are investigated using numerical examples and the results are compared with Monte Carlo simulation.Keywords Random eigenvalue problem · Asymptotic method · Stochastic dynamical systems · Probabilistic structural mechanics Nomenclature κ (r 1 ,r 2 ) jk (r 1 , r 2 )th order cumulant of jth and kth natural frequencies D (•) (x) Hessian matrix of (•) at x d (•) (x) gradient vector of (•) at x