## 4 Interval Equations System Solving  
In the paper, [PERSON]'s fixed-point theorem was applied to solve the system of linear interval equations ([PERSON] [47], [PERSON] [58]):  
\[\mathbf{D}^{I}(\lambda)\widetilde{\mathbf{q}}^{I}(\lambda)=\widetilde{ \mathbf{P}}^{I}(\lambda) \tag{20}\]  
after substituting the appropriate vector \(\widetilde{\mathbf{P}}^{I}(\lambda)\). Equation (20) can be transformed into a fixed-point equation in the following form:  
\[g(\widetilde{\mathbf{q}}^{I})=\widetilde{\mathbf{q}}^{I} \tag{21}\]  
where  
\[g(\widetilde{\mathbf{q}}^{I})=\widetilde{\mathbf{q}}^{I}-\mathbf{R}\left( \mathbf{D}^{I}\widetilde{\mathbf{q}}^{I}-\widetilde{\mathbf{P}}^{I}\right)= \widetilde{\mathbf{R}}\widetilde{\mathbf{P}}^{I}+\left(\mathbf{I}-\mathbf{R} \mathbf{D}^{I}\right)\widetilde{\mathbf{q}}^{I} \tag{22}\]  
and \(\mathbf{R}\) is a nonsingular matrix. In Equation (22), all quantities except \(\mathbf{R}\) are interval numbers. Using [PERSON]'s fixed-point theorem, we can express the following:  
\[\widetilde{\mathbf{R}}\widetilde{\mathbf{P}}^{I}+\left(\mathbf{I}-\mathbf{R} \mathbf{D}^{I}\right)\widetilde{\mathbf{q}}^{I}\in\widetilde{\mathbf{q}}^{I} \tag{23}\]  
The solution to Equation (23) also satisfies the original Equation (20). Therefore, it is possible to write the following (Rump [59]):  
\[\widetilde{\mathbf{R}}\widetilde{\mathbf{P}}^{I}+\left(\mathbf{I}-\mathbf{R} \mathbf{D}^{I}\right)\widetilde{\mathbf{q}}^{I}\subseteq int\left(\widetilde{ \mathbf{q}}^{I}\right) \tag{24}\]  
where \(int\left(\widetilde{\mathbf{q}}^{I}\right)\) denotes interior of \(\widetilde{\mathbf{q}}^{I}\). If we denote that \(\widetilde{\mathbf{q}}^{I}=\widetilde{\mathbf{q}}^{C}+\hat{\mathbf{q}}^{I}\), where \(\widetilde{\mathbf{q}}^{C}\) is the computed vector for the central values of the parameters, and \(\hat{\mathbf{q}}^{I}\) is the uncertainty of the solution, then Equation (24) can be expressed as follows (Neumaier [45]):  
\[\widetilde{\mathbf{R}}\widetilde{\mathbf{P}}^{I}-\mathbf{R}\mathbf{D}^{I} \widetilde{\mathbf{q}}^{C}+\left(\mathbf{I}-\mathbf{R}\mathbf{D}^{I}\right) \hat{\mathbf{q}}^{I}\subseteq int\left(\hat{\mathbf{q}}^{I}\right) \tag{25}\]The equation used for the iterative determination of the uncertainty of the solution can thus be written as follows:  
\[(\hat{\mathbf{q}}^{I})^{n+1}=\mathbf{R}\widetilde{\mathbf{P}}^{I}-\mathbf{R} \mathbf{D}^{I}\widetilde{\mathbf{q}}^{C}+(\mathbf{I}-\mathbf{R}\mathbf{D}^{I} )\,(\hat{\mathbf{q}}^{I})^{n} \tag{26}\]  
According to work [45] by [PERSON], the matrix \(\mathbf{R}\) is defined as \(\mathbf{R}=\left(\mathbf{D}^{C}\right)^{-1}=\mathbf{H}^{C}\). The stopping criterion for the iteration is as follows:  
\[(\hat{\mathbf{q}}^{I})^{n+1}\subseteq(\hat{\mathbf{q}}^{I})^{n} \tag{27}\]  
The final value of the solution is obtained based on the following equation:  
\[\widetilde{\mathbf{q}}^{I}=(\hat{\mathbf{q}}^{I})^{n+1}+\widetilde{\mathbf{q} }^{C} \tag{28}\]  
The above algorithm can be successfully applied to solve systems of linear interval equations, as demonstrated in [43] by [PERSON] and [PERSON]. 5 Application of the Presented Method to Multi-Degrees-of-Freedom Systems with Viscoelastic Elements  
The matrix-form equation of motion for systems with viscoelastic elements is given in the form (12). However, such representation leads to interval quantities associated with one physical quantity appearing more than once, which can cause significant overestimations. As demonstrated [PERSON] and [PERSON] in [43] for statics and [PERSON] et al. in [44] for dynamics, overestimation is reduced by applying the element-by-element technique. This approach is detailed in [44] by [PERSON] et al. As a result of employing this technique, we obtain a singular stiffness matrix. Hence, it is necessary to introduce an additional constraint matrix \(\mathbf{X}\) for the considered system and apply the Lagrange multiplier method \((\beta^{I})\). Then, the equation of motion looks as follows:  
\[