The paper presents a method for determining the dynamic response of systems containing viscoelastic damping elements with uncertain design parameters. A viscoelastic material is characterized using classical and fractional rheological models. The assumption is made that the lower and upper bounds of the uncertain parameters are known and represented as interval values, which are then subjected to interval arithmetic operations. The equations of motion are transformed into the frequency domain using Laplace transformation. To evaluate the uncertain dynamic response, the frequency response function is determined by transforming the equations of motion into a system of linear interval equations. Nevertheless, direct interval arithmetic often leads to significant overestimation. To address this issue, this paper employs the element-by-element technique along with a specific transformation to minimize redundancy. The system of interval equations obtained is solved iteratively using the fixed-point iteration method. As demonstrated in the examples, this method, which combines the iterative solving of interval equations with the proposed technique of equation formulation, enables a solution to be found rapidly and significantly reduces overestimation. Notably, this approach has been applied to systems containing viscoelastic elements for the first time. Additionally, the proposed notation accommodates both parallel and series configurations of damping elements and springs within rheological models.