## 8 Discussion  
The methodology presented in this paper enables the calculation of arbitrary-order sensitivities of eigenvalues and eigenvectors using HYPAD. In this case, multidual numbers were used to compute the sensitivities with high accuracy. The methodology is grounded on the works presented by [PERSON] and [PERSON] [17] and [PERSON] and [PERSON] [33] and formulated for self-adjoint systems. Although the focus of this work was on the analysis of structural dynamical systems, without loss of generality, the same methodology could be applied to other types of physics described by the GEP (e.g., Buckling analysis, principal component analysis (PCA), etc.).  
The method was developed for the case of distinct and repeated eigenvalues and showed high accuracy for both types of problems when compared to repeated form analytical solutions in free vibration analysis of a cantilever beam. In the mathematical development of the equation to compute eigenvectors' sensitivities [33], it is stated that the matrix \([\Theta]\) can become singular in some cases when eigenvalues of high multiplicity (repeated, with several similar eigenvalues) are obtained. Despite this, no sign of singularity has been observed in the verification analysis performed in this work. However, a more in-depth analysis is needed in future to present a more in-depth conclusion.  
The implementation of the current approach into commercial structural analysis software (i.e., Abaqus, Ansys, etc.) is left as future work. This could allow an increase in computational efficiency and enable the widespread use of the method for general structural systems regardless of the number of degrees of freedom. Although multidual numbers were used in this research to eliminate both truncation and subtractive cancellation errors, traditional complex numbers can be used if the user requires only the computation of first-order sensitivities, which are available in most of the commercial programming platforms. The only consideration needed is that the perturbation step size \(h\) should be small, say \(\approx 10^{-10}\) times the variable of interest.  
The proposed methodology is limited to analyzing the linear dynamics response of structural systems and does not consider the effect of material non-linearities, damping, loading or geometrical non-linearities. In addition, the computational cost, in terms of time and memory consumption, is acknowledged as a limitation of the proposed methodology. This limitation arises from the loss of sparsity and memory when evaluating the matrix of the coefficients in Equation (4). The issue becomes particularly challenging in systems with a large number of input design parameters, as it requires solving numerous systems of equations to assess all combinations of higher-order derivatives. To address this limitation, we are actively developing a solution that involves implementing a local residual formulation capable of computing sensitivities from expressions derived at the elemental level. This new formulation aims to preserve the sparsity of the matrix of coefficients. By adopting this strategy, we anticipate significant improvements in the efficiency and scalability of our methodology.  
Figure 10: Step size refinement for the first mode of vibration. The dashed line represents the results from HYPAD and the solid line the results from FD. (**a**) First-order derivatives. (**b**) Second-order derivatives. (**c**) Third-order derivatives.