A proper generalized decomposition based Padé approximant for stochastic frequency response analysis

This article presents a proper generalized decomposition (PGD) based Padé approximant for efficient analysis of the stochastic frequency response. Due to the high nonlinearity of the stochastic response with respect to the input uncertainties, the classical stochastic Galerkin (SG) method utilizing polynomial chaos exhibits slow convergence near the resonance. Furthermore, the dimension of the SG method is the product of deterministic and stochastic approximation spaces, and hence resolution over a banded frequency range is computationally expensive or even prohibitive. In this study, to tackle these problems, the PGD first generates the solution of stochastic frequency equations as a separated representation of deterministic and stochastic components. For the banded frequency range computations, the deterministic vectors are exploited as a reduced basis in conjunction with singular value decomposition. Subsequently, the Padé approximant is applied based on the PGD solution, and the stochastic frequency response is represented by a rational function. Through various numerical studies, it is demonstrated that the proposed framework improves not only the accuracy in the vicinity of resonance but also the computational efficiency, compared with the SG method.

simulation, perturbation, Neumann expansion, and polynomial chaos expansion have been applied to the stochastic frequency response analysis. Monte Carlo simulation (MCS) 2,3 is a sampling-based technique that only requires evaluation of the deterministic solver over a set of input samples. However, since a large number of samples should be drawn to achieve sufficient accuracy, the computational cost becomes prohibitive when the underlying problem is already intensive. The perturbation method 4,5 utilizes the sensitivity of modal properties (i.e., natural frequencies and mode shapes), and Neumann expansion 6 expresses the inverse of the dynamic stiffness matrix as a convergent series. Although both methods are simple and computationally efficient, they can only be applied to a low level of uncertainty, and the stochastic response by the Neumann expansion can diverge near the resonance. 7 Polynomial chaos expansion (PCE), also known as the spectral approach, explicitly represents the stochastic response utilizing a set of orthogonal polynomials with respect to the input probability measure. 8 The coefficients of the PC basis functions are determined by projecting the residual of the model output onto the approximated stochastic space, and can be carried out through intrusive or nonintrusive ways. The stochastic Galerkin (SG) 9 method is an intrusive approach that solves enlarged equations by modifying the deterministic FEM codes. In nonintrusive approaches, the coefficients are obtained through the numerical quadrature 10 or least-square 11 methods. Although PCE-based methods converge faster than MCS under certain regularity conditions and are not limited in terms of the uncertainty level, the computational cost grows dramatically as the number of random variables or polynomial degree increase.
In the field of uncertainty quantification, the SG method provides a general framework for the numerical simulation of stochastic governing equations. 12 However, when applied to the stochastic frequency response analysis, it has been reported that the accuracy of the first two statistical moments is significantly degraded near the resonance. 13,14 From the mode superposition perspective, the stochastic response can be expressed as a rational function of the modal properties, and it becomes highly nonlinear (nonregular) around the resonance as the denominator approaches zero. Accordingly, the classical PC basis is inefficient to approximate these response characteristics with slow convergence, and several studies have been conducted to tackle this issue. Noting that the first two statistical moments (mean and standard deviation) show spurious oscillations, Aitken's transformation 15 was applied to accelerate the convergence. Multielement approaches 16,17 have been investigated for frequency response analyses, where the global response is constructed by the local PC expansion of decomposed stochastic space. Rather than utilizing the PC basis directly, the SG-based Padé approximant 18 and sensitivity-based stabilization 19 methods additionally employ a denominator polynomial of the random variables to reflect the nonregular response nature.
In general, there are three difficulties in applying the SG method to the stochastic frequency response analysis: (1) the high nonlinearity of the stochastic response with respect to the input uncertainties; (2) the dimension of algebraic equations, which is the tensor product of the deterministic (spatial) and stochastic approximation spaces; and (3) the excessive computation of the frequency equations over the banded range, which increases linearly with the number of frequency evaluations. In this study, to overcome these difficulties, a proper generalized decomposition (PGD) based Padé-approximant is proposed. The main contribution of this article is that it suggests an approach that combines the advantages of the PGD 20,21 to alleviate the dimensionality of the SG method, and the Padé approximant 18,22 to reflect the nonlinear nature of the frequency response. In the present PGD context, the solution of the frequency equations is constructed by the separated representation of the deterministic vectors and stochastic functions, and their components are obtained by the updated progressive Galerkin scheme. To accelerate the computational efficiency over the banded frequency range, the deterministic vectors are reutilized along with singular value decomposition. Finally, the Padé approximant represents the stochastic response by a rational function, and its coefficients are computed based on the PGD solution. The rest of this article is organized as follows. In Section 2, we present the background of the spectral stochastic frequency response analysis. Section 3 describes the proposed PGD-based Padé approximant. In Section 4, the proposed framework is verified through numerical studies. Finally, we provide concluding remarks in Section 5.

Discretization of random field
If a random field is considered as the input uncertainty, it should be discretized prior to performing the stochastic frequency response analysis. Let us consider the spatially correlated random field a(x, ) defined on a domain x ∈ D with probability space (Θ, B, P), where ∈ Θ is a sample from the sampling space Θ, B is the σ-algebra, and P is the probability measure, respectively. The mean a 0 (x) and covariance function C a (x 1 , x 2 ) of the random field are defined by where ⟨⋅⟩ is the expectation operator. According to Mercer's theorem, 9 the covariance function can be decomposed as where k , k are the eigensolutions of the covariance function, which are obtained through the Fredholm integral equation of the second kind as follows: Based on the mean and eigensolution, the random field can be represented by truncated Karhunen-Loève (KL) expansion as 1 If the random field is modeled as a Gaussian, the standard random variables with zero mean and unit variance are used. For a non-Gaussian random field, mean-square convergent series utilizing polynomial chaos can be applied. 8 After the modeling step (4), the random field is equivalently formulated by the random variables defined on the finite-dimensional probability space (Θ (M) , B (M) , P (M) ). 23 In obtaining the eigensolution of Equation (3), the analytical solutions only exist for simple geometries with limited types of covariance functions. Therefore, to extend the applicability, the finite element method (FEM) is adopted in this study. Based on the discretized geometry, the eigenfunction is approximated by (x) = N(x) , where N(x) and are a shape function matrix and eigenvector. Substituting this relation into Equation (3) and applying the Galerkin projection onto the residual, the following symmetric eigenvalue problem is obtained: where matrices C and B are computed by Once Equation (5) is solved, the truncation number M is selected by examining the decay rate of the eigenvalues, and the random field is finally represented by Equation (4).

Stochastic Galerkin method
In the present study, a linear dynamical system in the frequency domain is considered. Suppose that the input uncertainties are modeled as mutually independent random variables = { 1 , … M } defined on the probability space, and let S = L 2 (Θ (M) , B (M) , P (M) ) be the associated stochastic space of all square-integrable random variables. 20 The stochastic frequency equations of N degrees of freedoms (DOFs) utilizing the FEM can be expressed as where and i = √ −1 are the excitation frequency and the imaginary unit; M( ), C( ), K( ), A( , ), and f( , ) are the mass, damping, stiffness, dynamic stiffness matrices, and force vector, respectively. The stochastic response u( , ) is sought in the tensor product of N-dimensional spatial and stochastic spaces, denoted by C N ⊗ S. The dynamic stiffness matrix and force vector are determined depending on the modeling of input uncertainties. In this article, Young's modulus is assumed as a homogeneous Gaussian random field, and the relevant implementation details will be discussed in Section 4.
When approximating the stochastic space, polynomial chaos expansion (PCE) utilizes polynomial functions of random variables. The polynomial chaos (PC) basis is constructed starting from a set of univariate polynomials. Let Ψ k ( k ) be a polynomial of degree k which satisfies the orthogonality condition, that is, where k k and P( k ) are the Kronecker delta and probability measure of k . The classical families of univariate polynomials are determined based on the Wiener-Askey scheme. 8 The multivariate polynomials for random variables are then built by the tensor product of the univariate polynomials as where = { 1 , … , M } is a set of multi-indices. The most common way to truncate the multi-indices set is to set the maximum polynomial degree to p. The truncated set of Ψ ( ) is then denoted by The number of polynomials P is the cardinality of A M,p , which is given by For notational simplicity, a one-to-one mapping between the set A M,p and index (1, … , P) is applied so that polynomials are ordered from the lowest to the highest degree (i.e., the degree of Ψ i ( ) ≥ the degree of Ψ j ( ) for indices ∀i > j). Denoting the finite-dimensional space spanned by the PC basis as S P = span({Ψ i } P i=1 ) ⊂ S, the solution of Equation (7) is then sought in the tensor-product space of the spatial finite elements and stochastic multivariate polynomial C N ⊗ S P as follows: whereû j ( ) ∈ C N is the coefficient vector of each PC basis. Substituting Equation (11) into Equation (7) yields In the stochastic Galerkin (SG) method, 9 the coefficientsû j ( ) are obtained by projecting the residual of Equation (12) onto the PC basis Equation (13) is assembled into the system of linear equation as follows: where the submatrixÂ ij and subvectorf i are computed by.
Equation (14) has a block sparsity pattern resulting from the orthogonality of the PC basis and dynamic stiffness matrix, and each submatrix inherits from the sparsity pattern of the dynamic stiffness matrix. 24 Once Equation (14) is solved, response statistics can be easily estimated by performing a Monte Carlo simulation (MCS) on approximation (11).
Regarding the SG method, the dimension of the complex-valued equations for each frequency is NP, which strongly depends on the DOFs and the number of PC basis functions. For the random dynamical system, the stochastic response around the resonance becomes a highly nonlinear function of random variables, and a fine approximation of the stochastic space is required to reflect this response characteristic. 16 In this context, the resolution of Equation (14) can easily become intractable when the underlying problem (DOFs of the FE model) is already intensive.

Proper generalized decomposition
In the following, a proper generalized decomposition (PGD) is presented to alleviate the dimensionality problem arising from the SG method. The PGD method, also known as a separated representation, constructs the stochastic response of Equation (7) by a low-rank m approximation of the deterministic (spatial) vectorŵ n ( ) ∈ C N and stochastic function λ n ( , ) ∈ S P as Here, the componentsŵ n ( ), λ n ( , ), and rank m are not known a priori and should be determined under the given target accuracy. The stochastic function is expanded by the PC basis as follows: Comparing the approximation form of the SG and PGD methods, which are expressed by Equations (11) and (16), the classical SG method seeks the coefficients of each PC basis in the full approximation space C N ⊗ S P . On the other hand, the PGD method exploits the low-dimensional subspace by the most relevant deterministic and stochastic components, and it is expected to converge faster with the rank m ≪ P. 21 Substituting Equation (16) into Equation (7) gives The objective is then to obtain two components by minimizing the residual of Equation (18) in some sense (i.e., Galerkin projection, minimal residual formulation), 16 and in this article, the updated progressive Galerkin scheme is considered. The progressive approach assumes that the decomposition ∑ m−1 n=1ŵ n ( )λ n ( , ) of rank (m − 1) is known from the previous iteration, and a new enrichment couple (ŵ m , λ m ) is found at the current step. Inserting this relation into Equation (18) yields where r is the residual vector, which is computed by Let us denote the tangent manifold of the couple (ŵ m , λ m ) by (ŵ * m λ m +ŵ m λ * m ), whereŵ * m ∈ C N and λ * m ∈ S P are the test function of the spatial vector and stochastic function belonging to the same approximation spaces ofŵ n and λ m . Imposing the Galerkin projection of Equation (19) onto this manifold, the following two equations are obtained 21 : where the superscript H represents the conjugate transpose. Equations (21) and (22) are so-called deterministic and stochastic problems, which are nonlinear due to the product of two componentsŵ m and λ m . The natural algorithm to resolve these equations in the PGD framework is the power iterations. It consists of solving the problems for each approximation space separately, assuming the other component is fixed. Once the power iterations are resolved, the update problem is additionally considered to accelerate the convergence. Hereinafter, each problem and its characteristics will be discussed.

Deterministic/stochastic problems
In the power iteration, we first consider the deterministic problem (21) based on the stochastic function λ k m by where the superscript denotes the iteration count. The deterministic problem has the same structure of standard (nonparametric) frequency equations discretized by the FEM, and its dimension N is equal to the DOFs of the model. After the deterministic vectorŵ k+1 m is obtained from Equation (23), the stochastic problem (22) for each PC basis {Ψ i } P i=1 can be expressed as Equation (24) can be assembled into a system of linear equations as follows: where W = (ŵ  (25) is P, which is related to the stochastic space approximated by the PC basis. Given the maximum number of iterations q max , the two problems (23) and (25) are solved recursively until the following convergence criteria are satisfied: where tol,ŵ and tol, and are the tolerance for the deterministic and stochastic components, respectively, and || ⋅ || is the L 2 -norm measure. This criterion can be seen as the couple (ŵ k+1 m ,̂k +1 m ) is reaching the fixed values, and only a few iterations are required to converge in practice. 23

Update problem
After the power iterations are resolved, the update problem is considered to improve the convergence. The stochastic response of Equation (16) for the rank m approximation can be rewritten as are the sets of deterministic vectors and stochastic functions. In the update problem, the stochastic functions are computed based on the deterministic vectors as follows: Equation (28) can be seen as the SG problem projected onto the reduced basisŴ. To express the update problem (28) in a matrix form, let us denote the set of stochastic functions by ( , Substituting this relation into Equation (28), the update problem for each PC basis In the same manner as Equation (25), Equation (29) can be assembled into the system of linear equations with dimension mP as follows: where and are obtained by =Ŵ H

AŴ
and =Ŵ H f. The enrichment process ends when the PGD approximation satisfies the desired accuracy, and the following criterion is considered in this article: where tol,g is the tolerance for global iteration. Given the maximum rank m max , the power-iterations and its update problem are applied until criterion (31) is met, and a flowchart of the overall PGD method is presented in Algorithm 1. Concerning the computational aspects, the intensive steps in the PGD are the resolution of the deterministic, stochastic, and update problems. The equation sizes of the deterministic and stochastic problems are the dimension of the two discretized spaces N and P, while the dimension of the SG method is the product of the two spaces, NP. Although these two problems must be solved recursively, it was already shown that their computational cost is much lower than that of solving the SG method when the DOFs for FEM or the stochastic dimension are large enough. 21 Lastly, the dimension of the update problem is mP, and grows linearly as the rank increases. However, since rank m is much lower than the deterministic dimension m ≪ N, there is still a computational advantage over the SG method.

3.1.3
Application to the banded frequency range At each excitation frequency, the aforementioned PGD allows efficient construction of the solution based on the deterministic vectors and stochastic functions. However, for the banded frequency range applications, Algorithm 1 should be carried out for every frequency, and its computational cost increases with the number of frequency evaluations. In this study, a strategy 16 reutilizing the deterministic vectors as a reduced basis is adopted to alleviate this issue. Assuming that frequency sweep proceeds in ascending order from the minimum frequency, m deterministic vectors are obtained at the first frequency by Algorithm 1. Rather than generating the couples from the empty sets (Ŵ = {}, = {}) for the next unresolved frequency r , the stochastic response is first approximated by Here, the stochastic functions ( r , ) are determined by the updating problem (30) based on the deterministic vectorŝ W. The accuracy of the approximation (32) is then evaluated by the criterion (31), and if the criterion is met, the PGD computation at the frequency r is resolved. Otherwise, the enrichment steps starting from the response (32) are employed until the PGD approximation satisfies the given accuracy. Supposing that m e enrichment steps are employed to meet the criterion (31), the deterministic vectors are augmented byŴ =Ŵ ∪ {ŵ m+i } m e i=1 and are exploited for the next frequency evaluation.
Although the required enrichment steps are greatly decreased by reusing the deterministic vector, one remaining issue during the frequency sweep is that the dimension of the update problem (30) is monotonically increasing with the enrichment steps. To prevent excessive computation of the update problem, this article extracts the dominant feature of the stochastic response by utilizing the singular value decomposition (SVD). When the rank is greater than the given threshold m th during the frequency sweep, the truncated SVD is applied onto the coefficient matrixW where m t = min(m, N, and ∈ R N×m t denotes the left, right singular vectors, and a block diagonal matrix of the singular values s 2 k , respectively. Assuming that singular values are sorted by descending order, the smallest number m r that satisfies the criterion (34) is selected.
where tol,s is the SVD tolerance. The deterministic and stochastic components based on criterion (34) are then reconstructed intoŴ ≃Ũ m r ∈ C N×m r and ≃ m rṼ T m r ∈ (S P ) m r , and the dimension of the update problem can be efficiently reduced without significant loss of accuracy. Given n f frequencies { k } n f k=1 with ascending order, the pseudo-code for computing the stochastic response for banded frequency range is presented in Algorithm 2.

PGD-based Padé approximant
A previous study 13 revealed that the statistical moments estimated by the PCE exhibit slow convergence around the resonance. Since the PGD also utilizes the PC basis for the stochastic functions, additional consideration in terms of accuracy should be required. Motivated by the fact that the frequency response is a rational function of modal properties, this article considers the Padé approximant for the scalar response u( , ) as follows: where n j and d j are the PCE coefficient of the numerator and denominator; P n and P d are the numbers of PC basis functions for the maximum polynomial degree of numerator p n and denominator p d , respectively. In the deterministic framework, the Padé approximant has been mainly applied to fast frequency sweeps. 25,26 The rational function in this context is a monomial of the frequency, and its coefficients are obtained based on the Taylor series solution. In the present stochastic context, the numerator and denominator are expanded by the PC basis of the random variables, and its coefficients are computed based on the PGD solution as Multiplying the denominator in Equation (36) and projecting the residual onto the PC basis {Ψ i } P t i=1 for the maximum polynomial degree p t yields Let us assume that the degree p t is greater than p n (i.e., p t > p n ) and the PC basis is ordered from the lowest degree. Denoting the numerator and denominator coefficients by n = {n 1 , … , n P n } T ∈ C P n and d = {d 1 , … , d P d } T ∈ C P d , Equation (37) can be written as Here, the relation ⟨Ψ i Ψ j ⟩ = 0 for indices i > P n is utilized. Each component G 1 ∈ R P n ×P n , G 2 ∈ C P n ×P d , and G 3 ∈ C (P t −P n )×P d is then computed by In solving Equation (38), denominator coefficients d are first considered. Assuming the overdetermined condition P t − P n > P d of the matrix G 3 , the linear equations G 3 d = 0 can be reformulated into the following minimization problem 22 : Equation (40) can be easily solved by applying the singular value decomposition to the matrix G 3 , where the solution d is the right singular vector corresponding to the smallest singular value. Once the denominator is computed, the numerator coefficients are directly obtained by n = G −1 1 G 2 d from Equation (38). Finally, based on approximation (35), the response statistics are estimated by performing the MCS.
Regarding the derivation of the Padé approximant presented in this article, the four PC basis for the PGD, numerator, denominator, and projection should be determined by setting each maximum polynomial degree (p, p n , p d , and p t ). Among the various combinations, this study considers the same degree for both the numerator and denominator, and their degrees are chosen by p n = p d = ⌈p∕2⌉ − 1, where ⌈⋅⌉ is a ceiling function. Based on Equation (10), the overdetermined condition P t − P n > P d of the matrix G 3 can be written as .
The projection degree p t is then determined as the smallest degree that satisfies Equation (41) and the assumption p t > p n . Hence, once the PGD degree p is selected, the remaining three polynomial degrees p n , p d , and p t are automatically determined.

NUMERICAL EXAMPLES
Prior to validating the proposed method through numerical studies, some implementation issues and the development environments are briefly overviewed. In the following examples, Young's modulus E(x, ) is modeled by the homogenous Gaussian random field. Utilizing truncated KL expansion of Equation (4), the random Young's modulus is represented by where E 0 is the mean value of Young's modulus; M is the truncation number by selecting the largest eigenvalues of the covariance function and the 90% cut-off criterion is applied in this article. When discretizing the given geometry, the same finite element meshes for both random field and stochastic frequency response problems (Equations (5) and (7)) are applied. Given relation (42), the element stiffness matrix is expressed by K e = K e 0 + ∑ M k=1 k K e k , where its mean K e 0 and random part K e k are computed by where D(E 0 ) and B are the constitutive and strain-displacement matrices, respectively. Hermite polynomials are considered as the PCE basis for Gaussian random variables. Concerning the PGD, the following control parameters are considered: maximum number of power iterations q max = 10, maximum rank m max = 100, tolerance for power iterations tol,ŵ = tol, = 5 ⋅ 10 −2 , global iteration tol,g = 10 −3 , and SVD tol,s = 10 −3 . The rank threshold m th will be specified in each example. When choosing the initial stochastic function λ 0 m in the power iterations, the coefficient of the zeroth-order PC basis is one and the rest are set to zero. In estimating the response statistics, the MCS with 10,000 samples is considered as a reference, and the same samples are applied to stochastic Galerkin (SG), PGD, and PGD-based Padé approximant (PGD-Padé). Lastly, all the algorithms previously discussed in this study are developed with in-house code, written in MATLAB 2017Ra.

Cantilever plate
The first example is a cantilever plate, as illustrated in Figure 1. A total of 450 nodes and 406 quadrilateral elements are used to discretize a given geometry. The material parameters considered in this example are as follows: mean Young's modulus E 0 = 2.1 × 10 11 N∕m 2 , Poisson ratio = 0.28, density = 7700 kg∕m 3 , and thickness t = 3 mm. The covariance function of Young's modulus is assumed to take the following exponential form: where E , l x , and l y are the standard deviation and the correlation length in the x-and y-axis. Under the parameters E = 0.1E 0 , l x = 0.5, and l y = 1, the eigenvalue problem (5) is conducted. The eigenvalue ratio and the first three eigenvector results are presented in Figure 2. From the results, three eigensolutions (M = 3) are required to simulate the given random field. A harmonic load F = 1 N is applied to the corner of the structure, and the end edge of the y-axis is subjected to a fixed boundary condition. In formulating the plate element, the assumed natural strain (ANS) method 27 is adopted to alleviate the transverse shear locking. The stochastic frequency equations (7) is expressed as where A 0 ( ) = [− 2 M 0 + i C 0 + K 0 ] and K k are the mean dynamic stiffness and the random stiffness matrices, respectively. The damping is modeled as proportional damping C 0 = M 0 + K 0 , and the parameters = 10 −1 , = 10 −5 are applied. The deterministic dimension is N = 1305, and a frequency range 0-200 Hz with an interval of 0.5 Hz is considered (n f = 401). Within this range, five deterministic natural frequencies are included, and the results are listed in Table 1 with the modal damping ratios. Figure 3 represents the mean and standard deviation results estimated by the SG method. The quantity of interest (QOI) is the transverse displacement at the node where the harmonic load is applied. The considered polynomial degrees p are 3 and 4, and the stochastic dimensions for each degree are P 3 = 20 and P 4 = 35, respectively. The obtained results indicate that the accuracy near the resonance is significantly degraded due to the spurious oscillations, while the good agreement with the MCS results is observed except for this range. In particular, within the range of 130-200 Hz, the oscillations corresponding to the fourth and fifth deterministic natural frequencies are superimposed, deteriorating the accuracy over the broad frequency range.
Before applying the proposed PGD for the given frequency range, we first check the convergence of the PGD at two distinct frequencies 80 and 146 Hz using Algorithm 1. The frequency of 80 Hz is outside of the resonance zone, and 146 Hz is close to the fifth deterministic natural frequency. The considered PGD degrees p are 7 and 9, and the corresponding stochastic dimensions for each degree are P 7 = 120 and P 9 = 220, respectively. Besides the updated Galerkin PGD, the minimal residual PGD 16 is also investigated, and the relative residual with respect to the rank is plotted in Figure 4. The relative residual for the PGD approximation is computed by where the norm in Equation (46) is defined in the space C N ⊗ S P (i.e., ||f|| 2 = ⟨f H f⟩). The prefixes M and G denote the minimal residual and Galerkin PGD, and the suffix U indicates the presence of the update step. From the results, the updated Galerkin method shows better convergence than the other methods. Although the minimal residual methods show monotonic convergence as already studied, 28 poor convergence is observed especially at 146 Hz. If the update step is not employed, the convergence becomes worse for both methods. The PGD response is obtained by Algorithm 2 with parameter m th = 30, and the mean and standard deviation results are presented in Figure 5. Since the stochastic function is approximated by the PC basis, the inaccurate oscillation phenomenon in the SG method is still observed near the resonance. Based on the PGD response, the Padé approximant (35) is then applied, and the first two statistical moments are given in Figure 6. The polynomial degrees p = 7 and 9 for the PGD lead to Padé approximant degrees (p n ∕p d ) of 3∕3 and 4∕4, respectively. From the results, it is clear that the first two statistical moments are well estimated over the entire range even though the deterministic natural frequencies are adjacent. Therefore, by comparing the approximation forms of the SG, PGD, and PGD-Padé methods (Equations (11), (16), and (35)), it can be inferred that denominator polynomials play an important role in reflecting the stochastic response characteristics around the resonance. The probability density functions (PDF) of the stochastic response at frequencies of 80 and 146 Hz are estimated, and the results are given in Figure 7. All PDFs are obtained through the MATLAB function ksdensity, which is based on the kernel density estimation. From the results, it can be seen that all methods match well the MCS results at 80 Hz. However, only the PGD-Padé can capture the response nature at the frequency of 146 Hz around the resonance zone, unlike other methods. To assess the accuracy of the first two statistical moments quantitatively, the following average relative errors for the N DOFs response are considered 25 : where (⋅) app i ( j ) denotes the mean or standard deviation of the ith DOF at jth frequency; the superscript app represents the SG, PGD, and PGD-Padé methods. The average relative errors (47) of the first two statistical moments are given in Table 2. The obtained results confirm that the errors of the proposed method are remarkably reduced in contrast to the SG and PGD methods. Finally, to provide information on computational efficiency, the CPU time for the MCS, SG, and PGD-Padé methods is examined. In this study, the system matrices are assembled in a sparse format, and a multifrontal method (UMFPACK) is applied to solve linear equations. 29 The dimensions of the linear equations for the SG method with degrees 3 to 4 are NP 3 = 26,100 and NP 4 = 45,675. A comparison of the total CPU time for the 401 frequency points is presented in Figure 8. To analyze the reuse effect of the deterministic vectors combined with the SVD, we additionally consider the following two cases: 1) without SVD (PGD-Padé w/o SVD) and 2) without reusing the deterministic vectors as a reduced basis (PGD-Padé w/o RB). Although our implementation is not fully optimized, it is clear that the CPU time for the SG method is more expensive than that for the PGD-Padé methods in any case. This implies that even if the PGD method requires iterative resolution for deterministic, stochastic, and update problems, it still has a computational advantage over the SG method. Figure 9 shows the rank results for each PGD method during the frequency sweep. If the SVD is not applied, the accumulated numbers of bases for each degree 7 and 9 are 68 and 71, respectively. The obtained results indicate that the advantage of reutilizing the deterministic vectors is lost for degree 9 because the dimension of the F I G U R E 9 Rank of the PGD approximation in the cantilever plate update problem increases monotonically. However, the proposed method controls the number of bases by extracting the dominant modes, and its computational efficiency can be accelerated compared with the PGD-Padé w/o RB as shown in Figure 8.

Rotor system
The second example is a rotor system consisting of a shaft with two flexible discs. The finite element model is illustrated in Figure 10, and its geometry parameters are listed in Table 3. A total of 3520 nodes and 2048 hexahedral elements are applied, and the material parameters are considered as follows: mean Young's modulus E 0 = 2.1 × 10 11 N∕m 2 , Poisson ratio = 0.28, and density = 7700 kg∕m 3 . In this example, Young's modulus of the shaft is modeled as a random field, and the covariance function is assumed as where E and l x are the standard deviation and the correlation length in the x-axis. The eigenvalue problem (5) is conducted under the parameters E = 0.1E 0 and l x = 0.6. The obtained eigensolution results are presented in Figure 11, and three eigensolutions (M = 3) are selected based on the 90% cut-off criterion.
The system is modeled to rotate around the x-axis ( = Ω{1, 0, 0}), and a fixed boundary condition is applied at x = 0. An unbalanced mass m e = 0.1 kg is attached to the node at position x m = {−0.512, 0.2, 0}, as shown in Figure 10. To improve the element performance, the hexahedral solid element is formulated by the enhanced assumed strain (EAS) method. 30 The equation of motion is derived with respect to the fixed-reference frame, 31 and is expressed by where G 0 and f 0 are the gyroscopic matrix and force vector due to unbalanced mass, respectively. The damping is modeled as proportional damping C 0 = M 0 + K 0 , and the parameters = 1, = 10 −5 are used. The stochastic frequency equations (49) can be derived as [ where A 0 (Ω) = [−Ω 2 M 0 + iΩC 0 + iΩ 2 G 0 + K 0 ] is the mean dynamic stiffness matrix. The total DOFs of this FEM model is N = 10,368, and the frequency response is examined under the range of 1-150 Hz with an interval of 0.5 Hz (n f = 299). The first two statistical moments by the SG method are presented in Figure 12. The quantity of interest is the z-axis displacement at the node where the unbalanced mass is located. The polynomials degrees are p 3 and 4, and the stochastic dimension for each degree is P 3 = 20 and P 4 = 35, respectively. As already examined in the previous example, the results obtained by the SG method exhibit inaccurate oscillations around the critical speeds (i.e., natural frequencies) regardless of the polynomial degrees. However, for the frequencies that are distant enough from all critical speeds, the relative errors are decreased as the polynomial degree increases. Figure 13 shows the relative residual of the Galerkin and minimal residual PGD methods at two distinct frequencies 25 and 50 Hz. The considered PGD degrees p in this example are 7 to 9, and the stochastic dimensions for each degree are P 7 = 120 and P 9 = 220. At frequency 25 Hz, which is sufficiently far from all critical speeds, the Galerkin PGD shows better performance than minimal residual PGD regardless of the update steps. However, the convergence of the Galerkin PGD without updating is severely degraded at frequency 50 Hz, which is close to the second critical speed. This implies the updating steps in Galerkin PGD are crucial for convergence and must be taken into account.
The mean and standard deviation results obtained by the PGD and PGD-Padé methods are given in Figures 14 and  15. Since the DOFs in this example are much larger than the stochastic dimension, it is more advantageous to solve the update problem than to perform power iterations. Thus, this example uses the parameter m th = 50 when applying Algorithm 2. The PGD results again confirm that the PC basis is not suitable for reflecting the response characteristics around the resonance. Compared with the SG method, increasing the polynomial degree does not guarantee accuracy. However, the simple posterior correction based on the Padé approximant alleviates this inaccuracy without any prior information on the response nature, and the first two statistical moments are estimated with high accuracy for all frequencies within the range. The average relative errors (47) of the three methods are listed in Table 4. The obtained results indicate that due to the increased accuracy around the resonant zone, the errors of the PGD-Padé method are reduced by more than two orders of magnitude compared with the other methods. Therefore, this suggests that the proposed method can also be applied to the rotor problems with the gyroscopic effect, which causes asymmetricity of the dynamic stiffness matrix.
The CPU time results for 299 frequency points are presented in Figure 16. Regarding the SG method, another difficulty besides the CPU time is that memory is seriously demanded due to the size of the linear equations for each degree NP 3 = 207,360 and NP 4 = 362,880, respectively. However, the PGD alleviates these issues by decoupling the tensor-product dimensions and can be more applicable to high-dimensional problems than the SG method within limited computer resources. The rank results for PGD approximations are given in Figure 17. Because of the large difference between the DOFs and stochastic dimension, the effect of reusing the deterministic vector is dominant is this example. In the PGD-Padé w/o SVD, totals of 97 and 95 bases are required for degrees 7 and 9, and the equation sizes of the update problem are up to 11,640 and 20,900, respectively. Accordingly, the cost of the update problem becomes excessive as the frequency sweep proceeds. However, the proposed method efficiently compresses the number of bases through the SVD and enables better computational performance than the other methods.

Plate coupled with acoustic cavity
The third example is a plate coupled with an acoustic cavity, as illustrated in Figure 18. The upper part of the cavity is coupled with a simply supported flexible plate, and the other parts are modeled as a rigid boundary condition. For the structural part, 3520 nodes and 441 quadrilateral elements are used. The acoustic field is modeled as water, and 4851 nodes and 4000 hexahedral elements are applied. The material parameters are considered as follows: mean Young's modulus E 0 = 7.2 × 10 10 N∕m 2 , Poisson ratio = 0.3, structural density s = 2700 kg∕m 3 , thickness t = 2 mm, sound speed of fluid c a = 1500 m∕s, and fluid density a = 1000 kg∕m 3 . In this example, Young's modulus of the plate is modeled as a random field, which follows the covariance function of Equation (44). Based on the parameters E = 0.1E 0 , l x = 0.55, and l y = 0.45, the eigenvalue problem (5) is conducted. The eigensolution results of the problem (5) are given in Figure 19, and four eigensolutions (M = 4) are selected to simulate the random field. A harmonic load F = 1 N is applied to the node at the position x f = {0.35∕4, 0.29∕4, 0.14}. The governing equation of the acoustic field is parameterized by the velocity potential , 32 and the stochastic frequency equations of the structure-acoustic system are expressed as where A 0,s ( ) = [− 2 M 0,s + i C 0,s + K 0,s ] and A 0,f ( ) = [− 2 M 0,f + K 0,f ] denote the mean dynamic stiffness of the structure and acoustic domain; C 0,sa and K k,s are coupling and the random stiffness matrices, respectively. The damping is modeled as proportional damping C 0,s = M 0,s + K 0,s for the structural part, and the parameters = 0.5, = 2 ⋅ 10 −5 are used. The dimensions of the structural and acoustic domains are N s = 1243 and N a = 4851, and the total DOF becomes N = N s + N a = 6094. The considered frequency range is 30-230 Hz with an interval of 0.5 Hz (n f = 401). Six deterministic natural frequencies are included in this range, and results along with the modal damping ratios are listed in Table 5.
For the PGD degree p = 9, Figure 20 represents the relative residual of the Galerkin and minimal residual PGD methods at two frequencies 30 and 70 Hz. The frequency 30 Hz is far from the resonance region, and 70 Hz is close to the second deterministic natural frequency. The obtained results again confirm that only the updated Galerkin method can be applied to the frequency equations without prior information about the resonance region. The mean and standard deviation are investigated by the SG, PGD, and PGD-Padé methods. Figure 21 shows the first two statistical moments of the transverse displacement where the load is applied, and the average relative errors for both displacement and velocity potential are given in Table 6. The polynomial degree for the SG method is p = 4, and its stochastic dimension is P 4 = 70. When obtaining the PGD results, the parameter m th = 30 is applied, and the PGD-Padé with degree 4/4 is computed based on the PGD response. From the SG and PGD results, it can be seen that the inaccurate oscillations occurring in the previous structural problem are still present in the coupled structural-acoustic system. In particular, the magnitude of the statistical response within the range of 180-225 Hz is relatively smooth due to damping, but large relative errors are still observed. However, the PGD-Padé method precisely estimates that the statistical moments and its relative errors for both displacement and velocity potential are greatly reduced compared with the two methods. Lastly, the CPU time and rank of the PGD approximation for 391 frequencies are given in Figure 22. The dimension of the SG method is NP 4 = 426,580. For the PGD-Padé w/o SVD, the accumulated number of the reduced basis is 95, and its cost is more demanding compared with the PGD-Padé w/o RB. The CPU time result for the proposed method is consistent with the previous examples, confirming better computational efficiency over the other methods.

CONCLUSIONS
In this article, the PGD-based Padé approximant for the stochastic frequency response analysis has been presented. The proposed method mainly consists of two steps: (1) First, the PGD by the updated progressive Galerkin scheme is applied to solve the stochastic governing equations in a separate form of deterministic and stochastic components; (2) Second, the Padé approximant is utilized based on the PGD solution, and the stochastic response is represented by the rational function of the PC basis.
To validate the performance, three examples of a rectangular plate, rotor system, and structure-acoustic system have been considered. The convergence results of the PGD confirm that the updated progressive Galerkin PGD can be applicable to the frequency equations, and the update step plays an important role in its convergence. By comparing the first two statistical moments, the proposed method accurately estimates the response characteristics over the examined frequency range, while the SG and PGD methods using the conventional PC basis exhibit inaccurate oscillations in the vicinity of the resonance. This is attributed to the denominator polynomial of the Padé approximant, which reflects the highly nonlinear nature of the frequency response with respect to the input uncertainties. In addition to accuracy, the computational cost also has been investigated. Since the PGD requires the resolution of a few decoupled deterministic, stochastic, and update problems, it allows significant computational savings compared with the SG and MCS methods. Among the PGD-based methods, the SVD can further reduce CPU time by recompressing the deterministic vectors.
In this study, the PCE basis for the Padé approximant is determined by setting the maximum polynomial degree. Since various approximations can be generated from the combination of the polynomials, further research on optimal basis selection should be conducted. Moreover, mitigation of the stochastic dimension caused by a large number of random variables also should be studied to generalize the proposed framework.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.