Padé approximants and the modal connection: Towards increased robustness for fast parametric sweeps

To increase the robustness of a Padé‐based approximation of parametric solutions to finite element problems, an a priori estimate of the poles is proposed. The resulting original approach is shown to allow for a straightforward, efficient, subsequent Padé‐based expansion of the solution vector components, overcoming some of the current convergence and robustness limitations. In particular, this enables for the intervals of approximation to be chosen a priori in direct connection with a given choice of Padé approximants. The choice of these approximants, as shown in the present work, is theoretically supported by the Montessus de Ballore theorem, concerning the convergence of a series of approximants with fixed denominator degrees. Key features and originality of the proposed approach are (1) a component‐wise expansion which allows to specifically target subsets of the solution field and (2) the a priori, simultaneous choice of the Padé approximants and their associated interval of convergence for an effective and more robust approximation. An academic acoustic case study, a structural‐acoustic application, and a larger acoustic problem are presented to demonstrate the potential of the approach proposed.


Background
Most engineering problems to be solved by numerical methods for engineering design purposes are parameter-dependent, often requiring multiple analyses to be conducted with increasingly complex models. To reduce the associated computational cost, it is desirable to be able to evaluate the solutions within the full range of allowable variation of these design parameters with computationally adapted methods (eg, adapted to the smoothness of the solution, its range, … ). Frequency sweep analyses are an example of such parameter-dependent problems, where several direct solutions have to be calculated. To get the solution in the full spectrum of interest thus becomes prohibitively demanding when several analyses are required, such as may be the case for optimization purposes or for problems where an assessment of the influence of uncertainties may be necessary.
In the context of frequency-dependent problems solved by the finite element (FE) method, several approaches have been developed to limit the impact of their computational burden. Widely used for several decades, modal methods have proved to be powerful from both an experimental and a numerical perspective, providing insight into the dynamic behaviour of the system of interest, as well as a basis suitable for a range of reduced-order models (ROMs). [1][2][3][4][5][6] These approaches aim at establishing a reduced set of equations to be solved at all the frequency steps of interest in the spectrum. Despite recent efforts made to extend such methods to more general types of damping, [7][8][9][10][11][12] these are however mostly limited to be efficient for either conservative problems or problems including simple models of damping.
More recently, alternative developments have focused not only on reducing the size of the parametric problem of interest but also to account for the general parameter dependence that may be involved to widen the range of application. Among the major developments in this area, the use of Padé approximants, 13 in various forms, has proven to be very effective and has received considerable attention for fast frequency sweeps.
Lanczos-based interpolatory methods, such as the Padé-via-Lanczos algorithm 14 or its subsequent development for more general partial fields and unsymmetric systems with the matrix-valued Padé-via-Lanczos approach, or Krylov subspace methods, [15][16][17] have been extended to second-order problems, being however limited to constant non-proportional damping models. Additionally, these extensions are made by a linearization of the second-order problem which involves doubling the dimension of the state-space vector, thus substantially increasing the memory requirements. To alleviate this drawback, second-order Arnoldi-based methods have been proposed, 18 allowing for a structure-preserved reduced model, however still limited to constant non-proportional damping. It is only recently that methods have been developed for reduced-order models able to effectively account for general nonpolynomial frequency dependence of the linear system of equations. [19][20][21][22][23][24][25][26][27] They are based on the explicit calculation of the solution vector and its derivatives at a restricted number of main frequencies in the spectrum of interest. In Beattie and Gugercin, 22 these explicit derivatives are used to span a subspace, suitable for a reduced frequency interval in the spectrum, on which the global matrices are projected. Alternatively, as done in the present work, these successive frequency-derivative vectors may be used for a component-wise solution expansion via Padé approximants. [19][20][21][23][24][25]28 With this approach, an efficient computation of the solution and its derivatives is performed at a very restricted number of frequencies in the spectrum of interest, which allows to construct Padé approximants for the expansion of the solution in the vicinity of these main frequencies. This alternative, further referred to as the component-wise Padé approximants, additionally offers the possibility to directly target the solution reconstruction for a specific degree-of-freedom (DOF) subset of interest. While the literature demonstrates the efficiency of this method in general, as has also been confirmed by the author in previous works, few attempts have been made to select the set of main frequencies, other than by heuristic methods. The choice of these frequencies, however, plays a very important role from both an accuracy and efficiency perspective. The idea of an adaptive frequency windowing is first introduced in Tuck-Lee and Pinsky, 16 in the form of an iterative convergence scheme associated with the Padé-via-Lanczos algorithm. For each frequency step, within 1 frequency subinterval of approximation, an incremental error comparison between 2 reduced-order systems of increasing sizes is applied until a convergence criterion is met. Further works have focused on such iterative approaches targeting projection-based interpolatory methods, eg, using residuals compared between 2 consecutive steps of increasing size of the reduced system. 26,29 In Rumpler et al, 24 the authors presented an adaptive frequency windowing for the component-wise class of Padé approximant methods. This approach consists in a combination of sequential, a priori estimates of the intervals of convergence of the Padé approximants, with a posteriori refinements using an error estimator. A common deficiency of the methods mentioned above is the potential appearance of gaps of unconverged approximation between consecutive intervals. In fact, the Padé approximants are limited in their range of convergence, mainly for numerical reasons, and such gaps of convergence may imply a need to add intervals for only small frequency ranges. This lack of robustness, obviously strongly dependent on the problem of interest, may hamper the efficiency and elegance of Padé-based approximations.

Current work
A first step towards a robust Padé-based approximation is proposed, here presented in the context of fast frequency sweeps. It is based on the component-wise type of Padé approximations discussed above. The choice of this method is mostly motivated by its ability to target a subset of DOFs at which the solution is to be calculated and to provide piecewise analytical expressions of this solution. Although the efficiency of this approach may be discussed against that of projection-based methods for a range of specific problems beyond the scope of this contribution, these unique features of the component-wise approach are of specific interest here. The proposed approach consists in predetermining the subintervals of convergence following an a priori step (or offline step) where the poles of the response are computed. In this sense, the Padé approximants may be linked to modal methods. Key to the proposed approach, such a Padé-based method is built on the poles in the range of interest only, without the typical truncation error issues of concern with modal methods. Furthermore, as previously stated, it is suitable for a wider range of applications. In the scope of this paper, the study will be limited to conservative problems where the link between modal methods and Padé approximants is straightforward. Further work is under investigation to ensure a reliable extension to problems with more general types of frequency dependence.
In Section 2, this paper focuses on recalling the current state-of-the-art regarding the use of component-wise Padé approximants for fast frequency sweeps, together with some observations motivating the approach by the author. The proposed method is then presented in Section 3, including some general guidelines to ensure good convergence properties. In Section 4, the component-wise Padé-based approach is tested for its robustness on several applications: A small acoustic problem consisting of a rigid acoustic cavity allows to illustrate some aspects of the methods and their convergence properties, then, a structural-acoustic coupled problem and a large acoustic application complete the presentation at larger scales, including a discussion on approximations over multiple intervals, robustness, and scalability of the method.

Numerical procedure for the calculation of the Padé coefficients
The starting point of the component-wise Padé univariate sweep, as described e.g. in Rumpler et al, 24 is given by a linear system of N equations (and N DOFs) in the following form: where x is the independent variable corresponding to the parametric problem of interest, eg, the angular frequency for the examples discussed in the present work. In an FE problem, Z(x), U(x), and F(x) respectively represent the system matrix of the discretized problem, the solution vector and the vector of externally applied loads. A component-wise expansion of the solution vector may be sought as Padé approximants in the form where the solution vector U(x 0 ), of which u(x 0 ) is a component, is assumed to be known after solving the system in Equation (1) for x = x 0 . P L (Δx) and Q M (Δx) are 2 truncated power series in the variable Δx = (x − x 0 ), to the orders L and M, respectively, and given by The expansion in Equation (2) being done for a component u(x) of U(x) around x = x 0 , the subsequent scalar expressions in this section also correspond to scalar components of the solution vector U(x). In previous works, it was shown 13 that the coefficients of the power series in Equation (3) may be determined from the coefficients of the Taylor series expansion at order (L + M) where These coefficients p k and q k are indeed solutions of the system of linear equations resulting from equating the Padé approximant in Equation (2) to the Taylor series expansion Equation (4), such that RUMPLER After developing Equation (6), a set of (L + M + 1) equations emerges from the identification of the coefficients of equal order in Δx, and normalizing the zero-order denominator coefficient, ie, q 0 = 1, where This system of linear equations may then be solved in 2 steps. First, the denominator coefficients q k may be solved for with a system consisting of the lower subset of equations in (7), The numerator coefficients may subsequently be determined in a second step by simple algebraic operations, The calculation of these Padé coefficients, being dependent on the Taylor coefficients a k in Equations (4) to (5), relies on the ability to efficiently calculate the L+M successive partial derivatives of the solution vector, at the reference point x = x 0 . This may be achieved via a recursive scheme, using a Leibniz formula resulting from the differentiation of Equation (1) with respect to x, at order k, in x 0 , where the zero-order derivatives correspond to the non-differentiated functions, and the binomial coefficients are given by The recursive expression for U (k) (x 0 ) follows from extracting the highest-order term from the summation in Equation (10), This implies that the successive derivatives of U with respect to x, required for the determination of the Padé approximations, can be efficiently calculated as the solution of a full-sized system of equations, with multiple right-hand sides.

Comments on the procedure
In agreement with the review of the theory in the previous section, an efficient approach to calculate the Padé-expansion around a solution vector U(x 0 ) (only 1 reference point at x 0 is considered in the present discussion) essentially consists of 5 steps: 1. Factorize the system matrix Z of Equation (1) at x 0 , and calculate the solution U(x 0 ) (coefficient a 0 ). 2. Solve for the L + M successive derivatives of U with respect to x, at x 0 , using a recursive multiple right-hand-side procedure, Equation (12) (coefficients a 1 · · · a L+M ).

For each component of interest in the solution vector:
3.1 Solve for the denominator coefficients via a small linear system of equations, Equation (8) (coefficients q 1 · · · q M , q 0 = 1). 3.2 Evaluate the numerator coefficients via simple algebraic operations as presented in Equation (9) (coefficients p 0 · · · p L ). 3.3 Evaluate the solution expansion around x 0 by evaluating Equation (2).
Note that while the first 2 steps have to be performed at a global scale, involving the full size of the original problem, the last 3 substeps may be limited to the DOFs of interest for the solution. A few remarks may be made in relation with the 5 main steps detailed above: • Step 2 involves a multiple right-hand-side problem with as many recursions as the number of numerator and denominator coefficients. In practice, there is a limit to the number of recursions that can be performed: The propagation of the small numerical error made at each recursion eventually affects the accuracy of the successive partial derivatives. • More importantly, the system matrix involved in the solutions of Equation (8), in substep 3.1, has the form of a Toeplitz matrix, which, although benefiting from very efficient algorithms, may become rapidly ill-conditioned with its increasing dimension. This implies that the solution of Equation (8) may become sensitive to small errors of the right-hand side, which consists of the Taylor coefficients, depending on the accuracy of the successive partial derivatives. • In addition to the previous points, the right-hand side of Equation (8), as well as the components of the system matrix, consists of the Taylor coefficients of highest order, ie, those most costly to calculate (latest stages of the recursive procedure in Step 2) and those with the most cumulated approximation error.
Even though the procedure presented in this section has proved to be very efficient for a wide variety of examples, 20,24,25 it appears that some of its steps may be improved in light of the points aforementioned. The following section discusses one such possibility, which additionally opens for the possibility of an a priori choice of Padé approximants, their associated range of convergence, and an improved scalability for the solution expansion of multiple solution vector components.

On the convergence of Padé approximants
The convergence of Padé approximants is a broad topic which has been given much attention by specialists for decades. 13,30 In the present work, a key result from the so-called Montessus de Ballore theorem, on the convergence of a series of Padé approximants, 31 is highlighted. This theorem applies to assumed meromorphic functions within a disk of sought convergence, which is typically the case for the parametric FE solutions of interest in this work. It guarantees, provided a sequence of Padé approximants whose poles correspond to simple poles of the approximated function within the considered disk, the uniform convergence of this sequence on this disk for increasing numerator degree and fixed denominator degree. This result implies that, given a set of simple poles of the solution function defining a disk of anticipated convergence, uniform convergence will be observed on this disk except at the poles location, for a denominator degree corresponding to the number of poles and a high enough numerator degree.
This lies at the basis of the proposed approach, detailed in the following section. The focus is to determine, a priori, both the Padé approximants and their anticipated domain of convergence, thus opening for the systematic approximation of FE solutions involving parametric sweeps over a wide range.

An improved procedure for robustness and efficiency
Following the observation made on the convergence of Padé approximants, their poles play a key role to ensure a reliable reconstruction of the solution. In the context of the FE problems of interest in this work, the poles of these Padé approximants should correspond to the zeros of the characteristic equation of the FE system matrix. Thus, starting with a step at a global scale (full-sized system of equations), by determining the zeros of this characteristic equation within subsets of the range of parametric sweep, would allow to address the observations highlighted at the end of Section 2.2: • The component-wise step consisting in solving for the denominator coefficients with Equation (8) would be replaced by an initial, offline, step at the global level; • The same poles would be ensured from one component to the other of the approximated solution vector, and thus, in light of the convergence properties recalled in the previous section, the same domain of convergence for all components would be enforced; • The number of Taylor coefficients required to be recursively calculated with Equation (12) may be reduced or reallocated to have numerator polynomials of higher degrees, thus possibly enhancing the accuracy over the domain of convergence of the approximants; • The information provided by the location of these poles may be used to decompose the full range of the parametric sweep into subdomains corresponding to 1 reference point of solution expansion and its associated interval of convergence.
In the scope of this paper, the system matrix Z(x) in Equation (1) is supposed to be a quadratic, real-valued function of x with real-valued eigenvalues, solution of a polynomial characteristic equation. Thus, Step 3.1 in Section 2.2 may be replaced by a step at the global scale, including the calculation of the M roots {x 01 , · · · ,x 0M } closest to x 0 , solution of the characteristic equation The steps of Section 2.2, for the Padé-based approximation of the solution around x 0 , may therefore be updated to the following 6-step procedure aiming at the same domain of convergence (ie, the same number of poles): 1. Solve for a set of roots {x 1 , · · · ,x n } , n ⩾ M, solution of the characteristic Equation (13), including the subset of the M closest roots to x 0 denoted {x 01 , · · · ,x 0M }. 2. Factorize the system matrix Z of Equation (1) at x 0 , and calculate the solution U(x 0 ) (coefficient a 0 for all DOFs). 3. Solve for the L ′ (as opposed to (L + M) in Step 2 of the initial procedure) successive derivatives of U with respect to x, at x 0 , using an iterative multiple right-hand-side procedure corresponding to Equation (12) (coefficients a 1 · · · a L ′ for all DOFs). 4. Once, and valid for all the components of interest in the solution vector, identify the M coefficients q 1 , · · · , q M (assuming simple roots at this stage) such that (coefficients q 1 · · · q M , q 0 = 1, common to all DOFs).

For each component of interest in the solution vector:
5.1 Evaluate the numerator coefficients via simple algebraic operations as presented in Equation (9) (coefficients p 0 · · · p L ′ ). 5.2 Evaluate the solution expansion around x 0 by evaluating Equation (2).
Arguably, this updated procedure may be seen as a component-independent methodology, as Steps 1 to 4 do not depend on the components of the solution vector, and Substeps 5.1 to 5.2 may be reconsidered from a global scale perspective. In fact, once the denominator coefficients are determined, the numerator coefficients are evaluated as a linear combination of the Taylor coefficients (see Equation (9)) involving the denominator coefficients. Thus, Substep 5.1 may be easily vectorized such that Steps 5.1 to 5.2 may be combined in a vectorized expression of Equation (2) for the complete solution vector in the vicinity of x 0 , where is a matrix of dimension N × (L ′ + 1), whose columns correspond to the Taylor coefficients for all the DOFs, determined after the (L ′ + 1) recursive solutions of the system of Equation (10), • Q L ′ is a square, upper triangular, Toeplitz matrix of order (L ′ + 1), such that Note that the inversion in Equation (15) consists, in practice, of a trivial scalar inverse associated with the denominator polynomial evaluation. This is enabled by the proposed approach where a common set of poles is enforced to all components of the solution vector. Furthermore, the link between a Padé approximation of the complete solution vector and an N -sized subset { } of its N DOFs is very straightforward with a localization matrix N such that where Û (x 0 + Δx) is the partial field, Padé-approximated, solution vector at the subset of DOFs { }. Thus, N is a rectangular boolean matrix of dimension N × N with a unique unit value in each column corresponding to the DOFs in { }. In practice, Equation (16) consists in selecting the lines of [A 0 · · · A L ′ ] corresponding to the DOFs of interest in { }.

Expansion of a solution over multiple intervals
When the parametric range of interest becomes too large to be approximated with 1 interval only (eg, encompassing a large number of poles or containing higher-multiplicity poles), the procedure described in the previous section may be very well suited for improving the robustness of decomposing this range into several intervals. Indeed, the first step in this procedure, where a sufficiently large number of poles may be determined-potentially all the poles in the range of interest-allows for an a priori decomposition in subintervals of Padé-based approximations. These intervals may typically be defined by a fixed number of poles or by pushing poles with a higher-order multiplicity to their bounds. In contrast with the original procedure in Section 2.2, this latter point prevents the procedure from being potentially hampered by the appearance of a higher-multiplicity pole within an interval of reconstruction, which would concentrate the approximation at a singularity where no convergence can be expected in the sense of the Montessus de Ballore theorem. Consecutive intervals may be chosen with 1 common pole defining the parametric bounds for each interval of anticipated convergence. Thus, considering 3 consecutive intervals centred on the interval associated with x 0 , these may write As assumed in this example, the number of poles per interval may be chosen as a constant or with little difference from one interval to the next. The full parametric range of interest may be decomposed a priori in this way.

APPLICATIONS
Two sets of examples are presented in this section. First, a small acoustic problem is used to illustrate some convergence differences between the 2 procedures presented in the previous sections. Then, more complex applications, including a medium-sized structural-acoustic problem and a larger acoustic test case, will allow to illustrate the potential of the approach proposed in this paper for both a partial field solution and a complete solution. In the following are the guidelines and criteria respected which have proved to be a good tradeoff between accuracy and efficiency in the tests considered by the author (this tradeoff being potentially specific to each problem), • Predetermined number of poles per interval: 4 to 5; • Interval bounds: 1 common pole for 2 subsequent intervals (as commented in the analyses, if accuracy at the junction between 2 intervals is critical, an overlap consisting of 2 common poles may be enforced to ensure continuity of the solution); • Expansion points: close to the centre of the interval, away from a pole, and if possible with a balanced number of poles on each side; • Numerator polynomial degree: twice that of the denominator, ie, 8 to 10; • Additional consideration: In the presence of a cluster of poles, try to separate them over 2 subsequent intervals.

Convergence on 1 interval: small acoustic cavity
In this section, some points of comparison between the original procedure in Section 2.2 and the updated one are analysed on a simple validation case consisting of a conservative acoustic problem. A rigid acoustic cavity of dimensions 0.1×0.15× 0.25 m 3 is excited at the corner (0, 0, 0) with a time-harmonic acoustic point source for consideration in the range f ∈ [500, 2250] Hz. The mesh of this small test case, suitable for the entire frequency range, consists of 4576 DOFs. The acoustic pressure fluctuation at the arbitrary position (0.06m, 0.11m, 0.16m) is used to compare the frequency sweeps obtained by using a Padé approximant on 1 interval with the reference solution. There are 11 eigenfrequencies of the acoustic cavity within the range of interest, see Table 1.
First, anticipating a limitation in the Padé approximants orders mostly due to the loss in precision associated with the recursive multiple-right-hand-side solution step, a convergence test is made for a fixed denominator degree of 7 (at most 7 poles in the solution can be captured). The 2 procedures are compared, for 1 Padé approximant with increasing orders of its numerator, which should imply a uniform convergence on an interval encompassing the 7 closest poles to the reference frequency. The reference frequency is first chosen at 1500 Hz, thus limiting the convergence to an interval including the eigenfrequencies {1143 Hz · · · 2058 Hz}. Figure 1 compares the convergence for these 2 procedures with Padé approximants of numerator expansion orders 6, 9, and 11. The main conceptual difference between the 2 procedures can be clearly seen in Figure 1A,B as the original procedure relies on the scalar-component data to determine the pole of the Padé approximant when the updated approach relies on the eigenfrequencies of the full-size problem. Consequently, the latter approach enforces the position of the poles no matter which step of convergence is reached; thanks to the polynomial order of the numerator. Thus, as early as expansion order L = 6 is reached for the numerator, the 7 expected poles are visible on the approximated solution of the updated procedure in Figure 1B. In contrast, it appears that 2 poles are missing for the original procedure (at 1847 and 2058 Hz). A typical behaviour is for the poles to be shifted away in such circumstances, in the higher-frequency region in this case. Upon convergence with increasing numerator order, one of the missing poles is shifted down to the domain of interest, as can be seen in Figure 1C,E. It however does not match the eigenfrequencies of the cavity, and 1 pole remains uncaptured by the approximation. The updated procedure, however, manages to capture the dynamic content of the response in the entire interval, all 7 eigenfrequencies as well as the sound pressure level in-between being accurately approximated after an expansion of the numerator to polynomial order 11.
The convergence limit of the original procedure is further illustrated in Figure 2, where 9 poles are considered, with a reference frequency of 1425 Hz, thus limiting the a priori convergence to an interval including the eigenfrequencies {686 Hz · · · 2060 Hz}. Note that this interval includes 2 very close eigenfrequencies (2058 and 2060 Hz) at 1 bound of the domain, a challenge for the approximation. Upon reaching limitations due to a combination of the high order of recursive derivatives in Equation (12) and the ill-conditioned nature of the problem in Equation (8), the Padé-approximation based upon the original procedure fails to capture the eigenfrequency at 686 Hz and collapses in the region 1715 to 1847 Hz and above (see Figure 2A,C). Indeed, estimating the location of the 9 poles captured by the approximant (L = 12; M = 9) in Figure 2C gives 2 complex conjugate poles together with real poles at {1143, 1210, 1333, 1374, 1717, 1809, 2431} Hz. Note that only 4 eigenfrequencies are accurately captured as poles of the approximant. In contrast, the updated procedure allows to have a good representation over the entire domain anticipated, even giving a fair approximation of the response at the challenging upper bound, where 2 eigenfrequencies almost coincide (see Figure 2D).

Multi-interval sweeps: structural-acoustic test case
In this section, a slightly larger, more complex problem is considered to test the proposed method on multiple intervals and to provide qualitative estimates of the performance. The model consists of a structural-acoustic example, including a rigid  acoustic cavity, except for 1 flexible wall. This structure is excited by a time-harmonic normal unitary force. The output is taken at a point within the acoustic cavity or as an average over all the points in the cavity (eg, mean quadratic pressure in the acoustic cavity). The details of the geometry, the mesh of the cavity involving 22 983 DOFs, and the arbitrary input and output points coordinates are provided in Figure 3. The frequency range of interest for this test case is limited to the range f ∈ [0, 200] Hz, which includes 46 eigenfrequencies. A reference solution of the acoustic level in this range, at the output point, is presented in Figure 4.
The position of the eigenfrequencies is also plotted. This reference solution shows not only the increased range of dynamic behaviour, which implies the need to approximate the solution over several intervals, but also the potential complexity to choose these intervals and ensure a good approximation. Two aspects in particular may be highlighted: (1) The modal density is very dependent on the location in the frequency range (eg, see the eigenfrequency gap around 150 Hz), which makes a fixed interval size for the approximation not an optimal choice, and (2) some frequency ranges locally have a high modal density with clusters of eigenfrequencies which may present a challenge for the Padé approximation. Due to these regions of high modal density, a quite fine frequency stepping is required, chosen here to be at intervals of 0.2 Hz.
The approximation, using the proposed method based on the use of Padé approximants, is calculated according to the guidelines at the end of Section 3.3. Following a preliminary step, involving the solution of the eigenvalue problem for   Table 2.
The necessary information for each interval consists in the main expansion frequency, around which the Padé expansion is calculated, together with the bounds for each interval. It is recalled that these bounds consist of common eigenfrequencies between 2 consecutive intervals. The Padé approximation corresponding to the procedure involving Steps 1 to 6 in Section 3.2 produces the results plotted in Figure 5.
The approximation shows, overall, an excellent match with the reference solution, even in challenging regions where clusters of eigenfrequencies are observed. A few local mismatches appear, eg, at 90.7 Hz (a structure-driven resonance)   and around 184 to 186 Hz, which are however located close or at the boundaries, and at eigenfrequencies, where there is no solution to the conservative problem, and no local convergence expected of the Padé approximants. A qualitative estimate of the efficiency of the method is done for an implementation in the same environment (Matlab, single threaded), using very close code architectures, and performed in a single run on the same machine. The performance may be summarized as follows: • Complete frequency sweep (1000 frequency steps): 25 minutes CPU time, ie, 1.47 seconds per frequency step.
• Padé approximation: approximately 6.6 seconds for the preliminary eigenvalue problem; approximately 34.5 seconds for the solutions and successive derivatives at the 12 main frequencies, ie, on average 2.9 seconds per main frequency; and approximately 0.05 seconds for the evaluation of the numerator Padé coefficients and the solution. Overall, this corresponds to approximately 41.15 seconds for the Padé-approximated solution.
To put these estimates into perspective, the overall time for the Padé-based approximation, is equivalent to performing less than 28 calculations of the full frequency sweep, ie, less steps than there are modes in the frequency range of interest.
Furthermore, following the comments made at the end of Section 3.2, the solution at all DOFs in the problem becomes only a matter of post-processing with the proposed method. An evaluation of this feature is made, considering the output quantity of interest is the mean quadratic pressure in the acoustic cavity, ie, where P 0 = 20 Pa is the reference sound pressure in the air and F the acoustic domain of interest. The procedure is tested in the same conditions as for the single DOF solution, only adding a post-processing step for both the complete and Padé-approximated solutions where the mean quadratic pressure in the acoustic cavity is calculated for each frequency increment. The results are plotted in Figure 6. These are consistent with the results for the single DOF, the agreement between the complete solution and its Padé-based approximation being very good. The additional computational cost of evaluating the solution at all DOFs for the calculation of the mean quadratic pressure in the cavity at all frequency steps is an extra 1.8 seconds of CPU time, to be added to the 41.15 seconds estimated for the single-DOF approximation. This confirms that the proposed approach has the side benefit of being equally interesting for partial field reconstructions or complete solution approximations.

Scalability: larger acoustic test case
To illustrate the scalability of the approach proposed, a larger calculation is performed on a simplified geometry of an industrial model of a light train.
The interior volume is modelled as presented in Figure 7. The mean quadratic pressure in response to a point acoustic source in the volume is calculated in the frequency range f ∈ [8,150] Hz. The mesh, consisting of 389 166 acoustic pressure fluctuation DOFs associated with quadratic tetrahedral elements, is not refined enough for the frequency range of interest but provides a critical test case for an evaluation of the Padé approach against the reference frequency sweep. Despite a simplification of some geometrical features of the original model, the mesh nevertheless requires local geometrical refinements, mostly due to narrow areas, which result in a rapidly increasing size of the problem to be solved.   The approximation using the mixed modal-Padé approach proposed follows the guidelines introduced at the end of Section 3.3 and the introduction of Section 4 regarding the a priori decomposition in frequency intervals and choice of the expansion frequencies.
The preliminary, offline, step in the procedure, involving the solution of an eigenvalue problem, highlights a decomposition on the basis of 57 eigenfrequencies in the range of interest, leading to a suggestion of 16 intervals for the solution expansion. These, together with the main frequencies around which the solution is to be approximated, are presented in Table 3.
The value of an a priori decomposition is manifest in this example, as the rapidly increasing modal density implies that the first 2 intervals cover 50% of the frequency range of interest. The associated reference solution, together with its approximation using the mixed modal-Padé approach are plotted in Figure 8. The interval bounds (dot-dashed lines), their associated main frequency of expansion (red circles), as well as the positions of the eigenfrequencies (blue crosses), are also plotted. The response shows a sharp contrast in the modal density between the lower half of the frequency range where the first few eigenfrequencies are separated by 5 to 10 Hz and the higher half of the frequency range where the separation between eigenfrequencies varies between 0.1 to 6 Hz with an average of 1.5 Hz. Such a small separation is challenging for the Padé-based reconstruction, which benefits in this case from the a priori choice of main frequencies positioned away from poles of the response. Due to this small separation, the step for the frequency sweep is set to 0.2 Hz, which may still be too coarse in certain regions of the spectrum, but provides a good tradeoff in relation with the low and high frequency densities. As shown in Figure 8, the Padé reconstruction is robust and provides an excellent approximation of the reference solution. As for the previous test case, local peaks of error are located at external boundaries of subintervals of Padé expansion, ie, at the outer eigenvalues of these intervals. However, no global convergence in the sense of the Montessus de Ballore theorem is expected at these points, as the response is not defined and its pointwise error sensitive to slight frequency shifts in the approximation. In cases where the degree of accuracy at the connection between 2 intervals may be critical, an eigenfrequency overlap between consecutive intervals may be considered.  To provide further qualitative insight into the relative efficiency and scalability of the approach, a similar evaluation to the structural-acoustic test case is done and commented. The performance may be summarized as follows: • Complete frequency sweep (711 frequency steps): 3h26 CPU time, ie, 17.4 seconds per frequency step; • Padé approximation (see also third column in Table 4): approximately 91.7 seconds for the preliminary eigenvalue problem; approximately 7 minutes and 25 seconds for the solutions and successive derivatives at the 16 main frequencies, ie, on average 27.8 seconds per main frequency; and approximately 32.1 seconds for the evaluation of the Padé approximation at all DOFs, the evaluation of the numerator Padé coefficients being in practice negligible (0.03 seconds). Overall, this corresponds to approximately 9 minutes and 29 seconds for the Padé-approximated solution at all DOFs.
These estimates may be further put into perspective and related to the structural-acoustic test case for performance and scalability purposes: • The overall CPU time for the Padé-based approximation is equivalent to performing around 32 calculations of the full frequency sweep (which comprises 57 eigenfrequencies), which corresponds to around 4.5% of the reference CPU time. The equivalent was of about 28 calculations of the full system for the structural-acoustic case (which comprises 46 eigenfrequencies) or around 2.9% of the reference CPU time. Anecdotally, the associated ratios of these equivalent numbers of full steps to the number of eigenfrequencies in the spectrum is of 56% in the train application and 60% in the structural-acoustic test case. Note that the apparent efficiency of the method applied to the train application is dependent on the size of the mesh: a finer reference mesh, better suited for the full frequency range, improves the apparent efficiency of the proposed approach, as further detailed below and in Table 4. • The preliminary (offline) step, consisting in determining the eigenfrequencies in the spectrum, is equivalent in CPU time to calculating around 5.27 frequency steps of the full problem. This is in line with the previous test case where the CPU time associated with the eigenfrequencies corresponded to 4.4 frequency steps of the full problem.
• In both applications, the CPU time associated with the preliminary eigenvalue problem corresponds to around 16% of the total Padé-based approximation time (16.1% for the train application and 15.4% for the structural-acoustic case). If the CPU time associated with the offline eigenvalue problem is equally distributed over all the subintervals of expansion, which is a fair assumption in the present examples, then the equivalent CPU time allocated to the calculation of the poles represents 20% only of the CPU time necessary to calculate the successive derivatives for each interval, for both examples. Therefore, the initial investment for the calculation of these poles, which is the originality of the approach, is well compensated by the points listed in Section 3.2, and in particular: reduction of the number of steps and number of successive derivatives necessary in the procedure, negligible time to calculate the Padé coefficients (0.03 seconds, overall, for the train application), increased robustness and a priori definition of the subintervals leading to a reduced number of these intervals.
A comparison with the original approach is also made for this train application, although only for qualitative purposes as part of the efficiency evaluation is dependent on the implementation. For the original Padé procedure, the 2 cases are considered where the solution is evaluated either at 1 DOF or for the complete problem. For the sake of comparison, the same intervals as those determined by the proposed approach are used, involving 57 poles and 16 intervals. The modal-Padé approach proposed is also tested on a refined mesh, suited for the frequency range of interest, and consisting of 2 031 961 DOFs. Table 4 summarizes these results.
Based on this estimated distribution of the CPU time, it may first be observed that the original Padé approach seems overall slightly more efficient than the modal-Padé approach proposed if the solution at a few DOFs is of interest (columns 2 and 3 in Table 4). However, it should be noted that this comparison is done on the basis of the a priori determined intervals of solution expansion, which are not available with the original method. Indeed, as the original method does not guarantee the convergence over the same poles for all DOFs, the actual solution calculated here was not converged for all intervals of solution approximation. In practice, as previously mentioned, the original method would require a heuristic approach for the choice of intervals, which for a given order of Padé expansion would increase the number of intervals, and thus the associated computation time. Additionally, the reasonable overall CPU time increase implied by the robust approach proposed (+22% to +29% for 1 DOF and all DOFs calculations, respectively, with reference to the original Padé for 1 DOF only) has the side benefit of a marginal increase of the CPU time if the solution at all DOFs is evaluated. In contrast, an evaluation of the solution for the full system with the original Padé approach more than triples the computation time. The proposed Padé approach for the full problem has an overall speedup factor of more than 21 for the 389 166-DOF test case. For the 2M-DOF problem, this overall speedup factor is around 24. The comparison with the original Padé method highlights the fact that the initial global pole localization step is well compensated by the subsequent robustness due to the a priori choice of intervals, and the versatility of the method equally suited for a partial field reconstruction or a full-sized problem solution evaluation.
This final application thus confirms the main objective achieved by the procedure proposed, bringing robustness to a Padé-based approach which, although it allows for piecewise analytical expansions of the approximated solution, may suffer from a lack of systematic procedure to choose the main frequencies of expansion and their associated intervals of convergence. Despite the relative subjectivity with which efficiency benchmarking may be performed (relation between mesh size and frequency content, choice in the order of Padé expansions, choice of Padé subintervals without a systematic procedure, degree of a priori knowledge of the problem, … ), the range of applications presented and the associated analyses support the fact that the efficiency and scalability of the Padé approach are maintained if not improved with the proposed procedure.

CONCLUSION
In this contribution, a Padé-based solution expansion method is proposed, providing increased robustness as compared with the current state of the art. The focus is specifically on the choice of Padé approximants, ie, the order of truncation of their denominator polynomials and, to a lesser degree, of their numerator polynomials, in connection with the choice of expansion points. It is shown, in particular, that an a priori decomposition in subintervals whose size is linked to the poles of the solution, follows the theoretical results of the Montessus de Ballore theorem on the uniform convergence of a series of Padé approximants. A 2-step procedure is subsequently proposed, consisting in a preliminary determination of the poles of the solution and the associated subintervals, followed by, for each subinterval, a Padé-based expansion of the solution. In this second step, the poles previously determined play a major role allowing to further ensure the efficiency and suitability of the method for both a partial field or a complete solution approximation. The examples presented allow to confirm the premises of the proposed method as well as its potential for a robust, efficient approximation of parameter-dependent solutions. Limited to conservative problems, these examples show that from the perspective of a comparison with modal methods for dynamic problems, the proposed approach only requires the eigenfrequencies included in the frequency range of interest, whereas modal methods usually need modes to at least twice the highest frequency of interest.
Two major areas are currently under investigation to further the development of the proposed approach. First, to increase the efficiency of the method, criteria for the truncation of the numerator polynomials and the refined positioning of the expansion points for each interval are being studied. These would allow to complete the guidelines provided in this contribution. Second, the extension to a wider range of parametric dependence is sought, eg, including non-polynomial parametric expressions and complex-valued solutions.