## 3 Sensitivity Analysis for Eigenfrequency Problems Using HYPAD  
The methodology used to obtain arbitrary-order sensitivities of the eigenvalues and eigenvectors in eigenfrequency problems using HYPAD with multidual numbers is based on a process with three stages, as shown in the schematic in Figure 1. The main inputsrequired to perform the sensitivities calculation process are the maximum order of the derivative \(m\) to be computed, and all the input parameters of the model \(=\{\alpha_{1},\ldots,\alpha_{j}\}\). In Stage 1, HYPAD is used to compute the system matrices \([K]\) and \([M]\), and their \(m^{th}\) order sensitivities with respect to any input model parameter. In Stage 2, the matrices \([K]\) and \([M]\) are used to solve the traditional real-valued GEP for structural systems. Finally, in Stage 3, the sensitivities of \([K]\) and \([M]\) obtained from Stage 1 (e.g., \(\left[\partial^{w}M/\left(\partial a_{1}^{b_{1}}\ldots\partial a_{j}^{b_{j}} \right)\right]\) and \(\left[\partial^{w}K/\left(\partial a_{1}^{b_{1}}\ldots\partial a_{j}^{b_{j}} \right)\right]\forall w=1,\ldots,m\)), and the eigenvalues and eigenvectors obtained from Stage 2 (e.g., \(\lambda_{i}\) and \(\phi_{i}\)), are fed into a recursive subroutine that uses HYPAD to systematically expand Equations (3) and (4) to enable the computation of all the sensitivities of the eigenvalues and eigenvector up to order \(m\). Below, a detailed description of each stage is presented:  
**Stage 1** is presented in Figure 2a and summarized as follows:  
1. The first step is to convert all the input parameters of the model into multidual variables of order \(m\). 2. Next, unitary perturbations, \(h\), are added to the independent imaginary axes of the multidual axis of the specific variables for which the sensitivities want to be calculated. Here, adding perturbations to the same parameter in two different independent imaginary axes results in having second-order derivatives with respect to the perturbed parameter. Similarly, perturbing two different variables in two different independent imaginary axes results in crossed second-order derivatives. Following this logic, any combination of arbitrary-order sensitivities can be obtained. 3. With the design parameters represented as multidual numbers and the appropriate perturbations applied, the structural system is discretized to obtain the global mass and stiffness matrices. Regardless of the discretization and spatial integration method used (e.g., FEM, BEM, SFEM, or others), when any input design parameter is represented as multidual and consequently perturbed, the mass and stiffness matrices will result in two \(n\times n\) multidual matrices \([M]^{*}\) and \([K]^{*}\). This approach is non-intrusive because the discretization and integration algorithms are unchanged; however, as the variables are multidual, algebraic operations are conducted using MultiZ. 4. The next step in the methodology is to extract the real and imaginary information of \([M]^{*}\) and \([K]^{*}\) to form \(p+1\) real matrices with dimensions \(n\times n\), with \(n\) being the number of DOF in the system, and \(p\) is the number of dual imaginary axes. Each matrix corresponds to an axis on \([M]^{*}\) and \([K]^{*}\). The matrices obtained with the information from the real axes of \([M]^{*}\) and \([K]^{*}\), correspond to the global \([M]\) and \([K]\) matrices for the system and are identical to those obtained from a traditional real-variable analysis; while the matrices formed with the imaginary axes (e.g., \([K^{e_{1}}]\),..., \([K^{e_{p}}]\), and \([M^{e_{1}}]\),..., \([M^{e_{p}}]\)) contain the information related to the partial derivatives of \([M]\) and \([K]\) with respect to the design input parameters that were perturbed. 5. Finally, to complete the first stage, the derivatives of the system matrices are computed following Equation (7). Figure 1: General methodology used to compute arbitrary-order sensitivities of the eigenvalues and eigenvectors using HYPAD.