Projection methods for stochastic dynamic systems: a frequency domain approach

A collection of hybrid projection approaches are proposed for approximating the response of stochastic partial differential equations which describe structural dynamic systems. In this study, an optimal basis for the approximation of the response of a stochastically parametrised structural dynamic system has been computed from its generalised eigenmodes. By applying appropriate approximations in conjunction with a reduced set of modal basis functions, a collection of hybrid projection methods are obtained. These methods have been further improved by the implementation of a sample based Galerkin error minimisation approach. In total six methods are presented and compared for numerical accuracy and computational efficiency. Expressions for the lower order statistical moments of the hybrid projection methods have been derived and discussed. The proposed methods have been implemented to solve two numerical examples: the bending of a Euler-Bernoulli cantilever beam and the bending of a Kirchhoff-Love plate where both structures have stochastic elastic parameters. The response and accuracy of the proposed methods are subsequently discussed and compared with the benchmark solution obtained using an expensive Monte Carlo method.


Introduction
The analysis of complex stochastically parametrised engineering structures has recently received significant interest. One of the main factors affecting the analysis is the computational cost associated with computing the response of a system. This can be mainly attributed to the dimension of the structure under consideration. In order address this issue, this paper proposes and compares a set of projection methods that approximates the response of dynamic structures with stochastic parameters.
In this work, stochastic linear damped structural dynamic systems are considered. The stochastic parameter associated with our governing hyperbolic partial differential equation can be characterised by a random parameter a(x, θ) on a bounded domain on R d and a probability space (Θ, F, P ), where θ ∈ Θ is a sample point from the sampling space Θ, F is the complete σ-algebra over the subsets of Θ and P is the probability measure. The displacement function of the system is given by u(x, t, θ). x ∈ R d represents a spatial position vector where d is the number of spatial dimensions and t ∈ R + represents the time. Through the use of well established stochastic finite element methods, a set of discretised linear equations can be obtained to represent the partial differential equations. Numerous methods have been suggested in order to solve or approximate the solution of the discretised set of equations.
Direct Monte Carlo simulation has been widely used in collaboration with stochastic finite element methods to generate the system response at an arbitrarily large number of sample points covering the input parameter space, following which the lower order statistical moments of the quantities of interest are calculated [1,2]. However this method is not favoured for the simulation of large systems. This is due to the convergence of the direct Monte Carlo simulation being slow with increasing dimension and associated variance of the input stochastic space. In addition to this, computing the exact solution of a high resolution finite element model at every sample point render the method extremely expensive. Several approaches have been proposed to reduce the computational effort. These include centroidal Voronoi tessellations [3], quasi Monte Carlo [4,5], Latin hypercube sampling [6], multilevel Monte Carlo [7], derivative-driven Monte Carlo estimators [8] and subset simulation [9].
Expansion based methods have been explored for approximating the response of the discretised set of equations. The perturbation method expands the stiffness and mass matrices, the forcing vector, and the response vector associated with the stochastic finite element method in terms of a truncated Taylor expansion [10,11]. However, the magnitude of the coefficient of variation associated with this method must be restricted [1]. If the coefficient of variation is increased, additional error may be induced. It has been shown that by combining a Taylor expansion with component mode synthesis, an efficient and accurate representation of the modes of a dynamic stochastic finite element structure can be obtained [12]. A method which is computationally competitive with respect to the perturbation method has been suggested by [13]. The proposed Approximate Principal Deformation Mode (APDM) approach utilises the fact that a systems' stiffness matrix can be expressed as a linear function of the uncertain parameters. Consequently, the APDM method can produce more accurate results for higher values of the coefficient of variation. This method has subsequently been used as a foundation for other uncertainty analysis methods [14,15]. The Neumann expansion method has also been applied to approximate the response vector [16]. This is achieved by expanding the inverse of the random matrix in a binomial type series. Numerous methods have been proposed based upon the Neumann expansion. [17] utilises the Neumann expansion to invert complex valued stiffness matrices whilst [18] proposes an acceleration technique in conjunction with the Neumann expansion. By utilising partial bivariate decomposition and successive matrix inversions [19], Ref. [20] has proposed a method to improve the computational efficiency of the Neumann expansion method.
Another class of methods for solving such problems are projection schemes. By utilising Wiener's 1938 work [21], [22] proposed a polynomial chaos expansion (PCE) for stochastic finite elements. This method produces a linear combination of Hermite polynomials and undetermined deterministic coefficients. A generalised form for the PCE was subsequently published by utilising functions from the Askey family of polynomials [23,24]. The PCE approach has been widely used in numerous fields including structural dynamics [25], heat transfer [26] and fluid dynamics [27].
Due to the large computational cost associated with the PCE, numerous cost saving approaches have been suggested [28,29]. Sparse polynomial chaos expansions have been wildly used as a surrogate for full models [30,31]. A stepwise regression method has recently been suggested to build a spare polynomial chaos [32]. The PCE approach has influenced many other recent works.
Ref. [33] has analysed the compressive sampling of polynomial chaos expansions, whilst a random discrete L 2 projection on polynomial spaces has also been proposed [34].
Comprehensive reviews of these techniques have been conducted in [40,41]. Ref. [42] and [43] have suggested adaptive reduced basis strategies in order to ensure a prescribed level of accuracy for a given output. These adaptive reduced basis methods have embedded a goal-oriented error assessment within their proposed methods. Techniques based upon stochastic response surfaces, where the quantities of interest are computed for certain θ ∈ Θ, have also received significant interest in recent years [44,45,46]. We refer the reader to [47,48,49] for a more comprehensive review of the field and the literature.
The random eigenfunction expansion method has been utilised in [17] to formulate a reduced random basis. In turn, Galerkin methods have been used in conjunction with eigenfunction projections to analyse the response of structures which are subjected to both a static load [50] and a dynamic load [51]. However questions still exist regarding the exact nature of these projections.
As a result, this study proposes a set of novel hybrid forward uncertainty propagation methods to perform a harmonic analysis of structural systems. Novel improved hybrid solution techniques have been proposed as an extension to the projection methods. These approaches have been incorporated by applying a multiplicative sample based Galerkin method. This addition aims to lower the induced error. All the proposed methods have been applied to two example problems i) the bending of a cantilever beam ii) the bending of a Kirchhoff-Love plate. Both structures have been discretised and include random parameters. The proposed methods are subsequently discussed and compared through the use of lower order statistical moments and appropriate relative error estimates. Through the use of the relative error estimates, the open questions that surround the hybrid uncertainty propagation methods can be addressed by identifying an optimal reduced projection.
A short overview of the stochastic finite element method is given in Section 2 and three projec-tion methods are presented in Section 3. Methods for reducing the computational effort associated with the proposed methods are discussed in the proceeding sections. Section 4 discusses a method for approximating the random eigensolutions associated with the methods, whilst Section 5 discusses a modal reduction. A sample based Galerkin error minimising technique is presented in Section 6 for each of the proposed methods. Expressions for the response moments of the different projection methods are discussed in Section 7. Section 8 provides a summary of the proposed methods. The approaches are then compared by applying the methods to a one dimensional stochastic Euler-Bernoulli cantilever beam and to a stochastic Kirchhoff-Love plate in Section 9. A summary of the results and the major findings are presented in Section 10.
2 The stochastic finite element method

Discretisation of the random field
In order to proceed with the stochastic finite element method, it is necessary to discretise the random field that is associated with the governing equation. We consider the stochastic parameter a(x, θ) to be a Gaussian random field with a covariance function C a : D × D → R defined on the domain D. The covariance function is positive definite, symmetric and square bounded. The random field a(x, θ) can be expressed by a truncated Karhunen-Loève series expansion. This expansion is achieved by performing a spectral decomposition of the covariance function of the field a(x, θ) = a 0 (x) where a 0 (x) corresponds to the mean function of the random field and ξ i are uncorrelated standard random variables. As the random field under consideration is Gaussian, the random variables are deemed as uncorrelated standard normal random variables with zero mean and a unit variance. λ i and φ i (x) correspond to the eigenvalues and eigenvectors satisfying the following Fredholm integral If the eigenvalues rapidly decay, the value of M could be kept relatively small in order to obtain an accurate depiction of the Gaussian random field. However as the correlation length of the process tends to zero, the number of terms required to obtain an accurate representation would increase.

Formulating dynamic systems in the frequency domain through finite element modelling
The methods of obtaining the discretised random form of the governing partial differential equations are well-established in stochastic finite element literature. By utilising the finite element method the equations of motion for a multiple-degrees-of-freedom structural vibration problem can be expressed as with the initial conditions set as In Equation (3) M(θ) and K(θ) denote the random mass and stiffness matrices respectively. C 0 and f 0 (t) denote the deterministic damping matrix and the deterministic applied force whilst t represents the time. The displacement is represented by u(t) and the first and second derivatives of the displacement with respect to time are represented byu(t) andü(t) respectively. The random mass and stiffness matrices can be expressed as follows In the above expressions, M 0 corresponds to the deterministic mass matrix and K 0 to the deterministic stiffness matrix. M i and K i are symmetric matrices which contribute towards the random components of M(θ) and K(θ). The random mass matrix has been modelled with p 1 random variables whilst the random stiffness matrix contains p 2 random variables. µ i (θ) represents the random variables associated with the random mass matrix, and ν i (θ) represents the random variables associated with the random stiffness matrix. ζ denotes a diagonal matrix which contains modal damping factors, thus It is assumed that the all the diagonal entries are equal, therefore ζ 1 = ζ 2 = · · · = ζ N . In order to satisfy this condition, the damping matrix C 0 takes the following form [52] We have made the assumption that the deterministic damping matrix belongs to the class of proportional damping matrices [52]. Therefore, the deterministic damping matrix has been simultaneously diagonalised with the mass and stiffness matrices by utilising the deterministic undamped eigenmodes. The detailed discussion of general non-proportional damping is beyond the scope of the present paper. In order to compute the dynamic response in the frequency domain, the Laplace transform of Equation (3) is considered. Taking the Fourier transform of Equation (3) results in Hereũ andf 0 are the dynamic response and the forcing in the frequency domain. The random variables associated with both the random mass and the stiffness matrices can be grouped so that ξ j (θ) = µ j (θ) for j = 1, 2, . . . p 1 and ξ j+p1 (θ) = ν j (θ) for j = 1, 2, . . . p 2 . In turn, Equation (9) can be re-written and expressed as where D 0 (ω) ∈ C N ×N represents the complex deterministic part of the system and D j (ω) ∈ R N ×N the random components. The total number of random variables, M , can be computed through summing p 1 and p 2 . For the given configuration, the expressions for D 0 and D j are as follows Therefore by combining the definitions of D 0 (ω) and D j (ω) with Equation (10), all the necessary components have been obtained in order to solve the discretised system of equations in the frequency domain. In the subsequent sections, different projection methods for efficiently approximating the response are compared. This is done for θ ∈ Θ and for every frequency value ω ∈ Ω. We refer the reader to [53] for a discussion regarding certain issues that incur during the application of the stochastic finite element method.

Derivation of the projection methods
We will initially consider three different projection methods. In order to compare the accuracy and effectiveness of the three proposed methods, a benchmark solution can be obtained by implementing for each frequency and realisation.
We aim to propose a set of methods which computes the response vector by projecting onto a basis with scalar coefficients. The rationale behind proposing different methods is to analyse the effect of the nature of the coefficients and their associated vectors. The first three methods under consideration have the following characteristics: • Projecting onto a stochastic basis with stochastic coefficients • Projecting onto a deterministic basis with stochastic coefficients • Projecting onto a deterministic basis with deterministic coefficients For all methods, we aim to keep the basis vector independent of the frequency. This is done in an attempt to reduce the computational effort if more than one frequency value were to be analysed.
We initially consider the case which incorporates the whole stochastic nature of Equation (10). For this method, we aim to represent the response by projecting onto a stochastic basis with stochastic where α j (ω, θ) ∈ C denotes the random scalars which are contained in α(ω, θ) ∈ C N , and a j (θ) ∈ C N denotes the stochastic basis. These basis are contained within the matrix a(θ) ∈ C N ×N . The values of α j (ω, θ) and a j (θ) can be obtained through numerous approaches. One such approach is by solving the following multi-objective optimisation problem While the above approach gives the generic framework for the evaluation of the α(ω, θ) and a(θ), the process can be computationally expensive due to the method's slow convergence rate. Furthermore the method could be numerically unstable as the solution may not be unique. In order to avoid calculatingũ DM CS (ω, θ), an expression for the above L 2 relative error can be obtained by observing the residual and by noting that the approximate error of the solution obtained when using Equation (14) isε Here the error measure is defined by using the DM CS approach as a benchmark solution. A closed form of the error in the domain space of D(ω, θ) can be obtained. The residual can be re-written as whereũ * (ω, θ) is the true solution of the system which cannot be evaluated exactly. We can treat the solution of the DM CS approach,ũ DM CS (ω, θ), as the benchmark solution. It is assumed the DM CS approach gives a better approximation of the true solution compared toũ 1 (ω, θ). Using e(ω, θ) =ũ 1 (ω, θ) −ũ * (ω, θ) as the true error, we write following Equation (18) D(ω, θ)e(ω, θ) = r(ω, θ) Thus the resulting true error vector is obtained as But e(ω, θ) cannot be computed exactly and we have to resort to the approximate error indicator.
We can define a bilinear form asD(a, b) = D(ω, θ)a(ω, θ), b(ω, θ) where ·, · denotes an inner product in L 2 (Θ) × R N . Hence, from Equation (20) we can deducē Using Cauchy-Schwarz inequality, we have where || · || E denotes the norm consistent with the bilinear formD(·, ·) on L 2 (Θ) × R N (analogous to the elastic potential energy norm for structural dynamic systems). Combining Equations (21) and (22) we obtain which indicates a lower bound for the true error e(ω, θ) in terms of the approximate error indicator ε(ω, θ). The equality holds only under special circumstances which have been detailed in [54].
However as the computation capacity required to implement such an approach is vastly higher than that required for the benchmark solution, a different approach is needed.

Projecting onto a stochastic basis with stochastic coefficients (M 1)
In order to implement a different approach, the generalised eigenvalue problem for the undamped case is initially considered where λ k (θ) and φ k (θ) are the kth undamped random eigenvalue and eigenvector. For convenience, matrices that contain the whole set of random eigenvalues and eigenvectors are defined as follows The eigenvalues are arranged in ascending order so λ 1 (θ) < λ 2 (θ) < . . . < λ n (θ) and their corresponding eigenvectors are mass normalised and arranged in the same order. It is apparent that As the undamped eigenvectors from a complete basis, it is possible to obtain the response of Equation (10) through projecting on the undamped eigenvectors. This can be done through using the above identities. The Laplace transform of Equation (3) is initially reconsidered The modal damping matrix is defined as follows where ζ corresponds to the diagonal modal damping matrix introduced in Equation (7). By using the following modal transformationũ(ω, θ) = Φ(θ)ȳ(ω, θ) and by pre-multiplying Equation (27) with Φ T (θ), we obtain By combining the modal damping matrix and the orthogonality relationships defined above, it can be shown that Then by inverting −ω 2 I + 2iωζΩ(θ) + Ω 2 (θ) , one has As −ω 2 I + 2iωζΩ(θ) + Ω 2 (θ) is a diagonal matrix, its inverse is easy to compute and computationally inexpensive. By pre-multiplying both sides of the above equation with Φ(θ), we have By reintroducingũ(ω, θ) for Φ(θ)ȳ(ω, θ) a dynamic response in the frequency domain can be This expression can then be rewritten as a summation, where N corresponds to the number of degrees of freedom associated with the dynamic structurẽ The response of the dynamic stochastic system under consideration has been represented in the same form as Equation (14). The random scalars, α j (ω, θ), correspond to the result of . In turn, these random scalars are projected onto the space spanned by φ j (θ).

Projecting onto a deterministic basis with stochastic coefficients (M 2)
Thus far one projection method has been proposed through Equation (34) where both the basis and the projection coefficients are stochastic in nature. If the value of N is large, the computational effort associated with computing the undamped random eigenvalues and eigenvectors can be considered very high. This is especially true if we sample for every θ ∈ Θ. In an attempt to lower the computational effort, we will consider a method that projects random scalars onto a deterministic basisũ The polynomial chaos approach is a method which projects onto a deterministic basis with stochastic coefficientsũ where H k (ξ(θ)) represents the polynomial chaoses (corresponding to the random scalars) and u k represents unknown deterministic vectors that need to be determined. The value of P is governed by the value of M and by the order of the polynomial chaos expansion. However the polynomial chaos approach is a notoriously costly method if the value of P or N is high. More importantly the basis u k is a function of ω, thus this method does not comply with the desired form stated earlier in this section.
Although mathematically erroneous, this paper proposes a method that combines undamped random eigenvalues with undamped deterministic eigenvectors. By exchanging the undamped random eigenvectors seen in Equation (34) for their deterministic counterparts, the response vector for this new method can be expressed as where φ 0j denotes the jth deterministic undamped eigenvector. Therefore we aim to see if the vast majority of the stochastic nature of the system can be incorporated by only using the undamped random eigenvalues.

Projecting onto a deterministic basis with deterministic coefficients (M 3)
As a worst-case scenario we consider the case of when all the eigensolutions are deemed deterministic. For this case, the response vector takes the following form where γ 0j (ω) ∈ C and c j ∈ C N are deterministic scalars and basis respectively. If both the undamped random eigenvalues and eigenvectors seen in Equation (34) are exchanged for their deterministic counterpart, the response vector can be expressed as where λ 0j and φ 0j denote the jth undamped deterministic eigenvalue and eigenvector respectively.
Due to all the terms in Equation (39) being deterministic, the stochastic nature of the system is not at all incorporated into the response vector. It can be deduced that Equation (39) provides the deterministic solution and therefore can be established as a worst case scenario. However if the coefficient of variation associated with the stochastic process is low, this method could provide an adequate approximation of the mean of the true solution.
If the matrices U DM CS ∈ C N ×m and U 3 ∈ C N ×m contain the solution vectors for all realisations of a given frequency for the benchmark and M 3 methods, the Frobenius norm of the relative error is given by where m corresponds to the number of realisations. If the matrices U 1 ∈ C N ×m and U 2 ∈ C N ×m are similarly defined for the M 1 and M 2 methods the following propositions can be made Due to both the eigenvalues and the eigenvectors retaining their stochastic properties in M 1 method, it is intuitively expected that the M 1 method would induce the least amount of error.
As the entire stochasticity of the response vector is expected to be captured by the stochastic eigenvalues in the M 2 method, this method is not expected to outperform the M 1 method. In a similar manner, as the M 3 is deemed to be a worst-case scenario, naturally this method will not outperform both the M 1 and M 2 methods.
At present, the computational time associated with both the M 1 and M 2 methods can be considered quite high, especially for a high degree of freedom finite element system. This is due to two reasons. The first being the large number of terms in the summations seen in Equations (34) and (37). At present, the number of terms in the series corresponds to the number of degrees of freedom. Secondly, calculating the random eigensolutions is computationally expensive. Combining these reasons with the need to simulate the methods for each θ ∈ Θ accumulates to a high computational effort. These issues have been addressed in the following section where approximation and truncation techniques are introduced.

Approximating the undamped eigensolutions
Calculating the exact undamped random eigensolutions can be extremely expensive, especially if the number of degrees of freedom is large. Thus a sensitivity approach to approximate the eigensolutions could computationally be a better option. The set of exact undamped random eigensolutions can be obtained through combining direct Monte Carlo simulations with the eigenvalue problem of the undamped system expressed by Equation (24). The set of random eigenvalues could consequently be used to obtain the natural frequencies of the system where ω k (θ) represents the kth random natural frequency of a realisation. Due to the low computational effort associated with the first order perturbation method, this method has been used to approximate both the undamped random eigenvalues and eigenvectors. The random eigenvalues can be approximated by the following equation where λ j0 is the jth deterministic undamped eigenvalue and dξ k (θ) a set of Gaussian random variables with mean zero and unit variance. The derivative of the undamped random eigenvalues with respect to ξ k can be obtained through differentiating and manipulating the eigenvalue equation denoted by Equation (24) [55]. This results in the following equation where λ 0j and φ 0j correspond to the deterministic undamped eigenvalues and eigenvectors. As the deterministic undamped eigenvectors are mass normalised, the denominator in the above equation In the instance of Equation (45), the values of both ∂M ∂ξ k and ∂K ∂ξ k are as follows where M k and K k − p 1 correspond to the random components of M(θ) and K(θ) introduced through Equations (5) and (6).
In a similar manner, the random undamped eigenvectors can also be expressed by a first-order where φ j0 is the jth deterministic undamped eigenvector and dξ k (θ) a set of Gaussian random variables with mean zero and unit variance. It is possible to obtain the derivative of the jth undamped random eigenvector with respect to ξ k by expressing the result as a linear combination of deterministic eigenvectors. This can be illustrated by The full algebraic detail of obtaining the derivative can be found in [55]. The final expression for ∂φ j ∂ξ k is given by The values of both ∂M ∂ξ k and ∂K ∂ξ k are identical to those given in Equation (46). This method requires all the deterministic eigenvalues and eigenvectors to be known. Furthermore the eigenvalues are required to be unique. For the case of repeated eigenvalues the proposed methods would still be valid, however a different method would be required to approximate the eigenvectors [56].
Thus by approximating the random eigensolutions it can be categorically concluded that the computational effort associated with each of the proposed methods are as follows where C represents the computational effort. Providing that the proposed methods for approximating the eigensolutions are accurate, by comparing Equations (41) and (50) it is apparent that a trade-off between the error and the computational effort is present. M 1 provides the most accurate representation of the response vector however its computational effort is the highest. M 3 provides the least accurate response, however this is achieved with a considerably lower computational effort.

Modal basis reduction
At present all the methods described in Section 3 require the calculation and summation of N terms.
In cooperation with the approximations seen in Section 4, these methods can be deemed rather inexpensive in comparison with the polynomial chaos approach. The polynomial chaos approach requires a set of N P algebraic equations to be solved where P corresponds to the number of polynomial chaoses. However, we aim to lower the computational effort of our proposed methods even further.
Through revisiting the ordering of the eigenvalues seen in Equation (25) it can be deduced that where λ j corresponds to the jth eigenvalue. From the scalar terms α j (ω, θ), β j (ω, θ) and γ 0j (ω) seen Equations (34), (37) and (39), it can be observed that the eigenvalues appear in the denominator.
The scalar α j is shown for illustration For the values of j satisfying λ j (θ) + 2i λ j (θ)ωζ > ω 2 , it is apparent that the value of the denominator increases as the value of j increases. The value of the numerator depends on the deterministic forcef 0 and the undamped eigenvectors. The numerator cannot be ordered in terms of magnitude for different values of j, however due to the mass normalisation of the undamped eigenvectors, it can be deduced that the value of the numerator does not vary significantly. Therefore it is established that the value of α j (ω, θ) generally decreases as the value of j increases. Consequently the upper limits of the summations seen in Equations (34), (37) and (39) can be lowered. In turn, these equations can be expressed as respectively, where n r < N ≪ N P . The value of n r can be defined in two ways (a) the value can be predefined depending on the system under consideration (b) by selecting a value for ǫ which is sufficiently small, n r can be selected such that λ 0 (nr ) is the largest deterministic eigenvalue that satisfies λ0 1 λ0 (nr ) > ǫ. If the accuracy of the truncated series is not sufficient, the accuracy can be improved by increasing the predefined value of n r or selecting a lower value for ǫ.

Sample based Galerkin error minimisation
Three different projection methods have been proposed in Section 3. The first projects random scalars onto a stochastic basis whilst the second projects random scalars onto a deterministic basis.
The third method projects deterministic scalars onto a deterministic basis. We have shown that it is possible to approximate the random eigensolutions that arise in the proposed methods in order to lower the computational effort. However these approximations, in addition to the modal reduction introduced in Section 5 introduces error into the calculation. This has motivated an error minimisation technique through applying a sample based Galerkin approach. As a result, in addition to the three projection methods introduced in in Section 3, the following three projection methods are proposed: • Galerkin approach with projecting onto a stochastic basis with stochastic coefficients (M 1G) • Galerkin approach with projecting onto a stochastic basis with deterministic coefficients (M 2G) • Galerkin approach with projecting onto a deterministic basis with deterministic coefficients The case of incorporating a sample based Galerkin approach with projecting onto a stochastic basis with stochastic coefficients is initially considered.

Galerkin approach with projecting onto a stochastic basis with stochastic coefficients (M 1G)
The response vector for the given case is modified to take the following series representatioñ Here α j (ω, θ) and φ j (θ) correspond to the random scalars and random eigenvectors seen in Equation (34) whilst c j (ω, θ) ∈ C are constants which need to be obtained for each realisation. This can be done by applying a sample based Galerkin approach. We initially consider the following residual where ξ 0 = 1 is used in order to simplify the summation. D i (ω), ξ i (θ) andf 0 (ω) correspond to the terms arising in Equations (9) and (10). By making the residual orthogonal to a basis function, the unknown c j (ω, θ) can be computed. As Equation (56) can be viewed as a projection onto a stochastic basis, the residual is made orthogonal to the undamped random eigenvectors As a sample based Galerkin approach is considered, applying the orthogonality condition results Through manipulating Equation (59) it is possible to re-write the equation in the following form By defining the vector c 1 (ω, θ) = [c 1 (ω, θ) c 2 (ω, θ) . . . c nr (ω, θ)] T , Equation (60) can be re-written as The number of equations that need to be solved in order to calculate the unknown vector c(ω, θ) corresponds to the value of n r . By increasing the number of terms from n r to n r + 1, the number of terms in Z 1 (ω, θ) increases by 2n + 1. Therefore the lower the dimension of the reduced system, the fewer the number of equations that need to be solved. This is of importance as the given procedure needs to be repeated for every realisation and for every frequency under consideration.

Galerkin approach with projecting onto a stochastic basis with deterministic coefficients (M 2G)
A similar approach can be implemented for the case containing undamped random eigenvalues and deterministic eigenvectors. The response vector for this given approach takes the following form Here β j (ω, θ) corresponds to the scalars introduced in Equation (37) and φ 0j to the deterministic eigenvectors also introduced in Equation (37). c j (ω, θ) ∈ C are the unknown constants that need to obtained for each realisation of each frequency. A similar sample based Galerkin method can be implemented to compute the unknown constants. In order to mimic the projection seen in Equation (37), the residual is projected onto deterministic eigenvectors rather than the random eigenvectors seen in Section 6.1. By following a similar approach to that seen in Section 6.1, the unknown constants c j (ω, θ) can be computed by solving the following set of linear equations.
Z 2 (ω, θ)c 2 (ω, θ) = y 2 (ω) j, k = 1, 2, . . . , n r where . . n r and c 2 (ω, θ) is a vector that contains the unknown constants c j (ω, θ) The number of equations that need to be solved to compute the unknown coefficients corresponds to the number of modes retained in the reduced models in Section 5. Similarly to the the method described in Section 6.1, this procedure needs to be repeated for every realisation and for every ω ∈ Ω.

Galerkin approach with projecting onto a deterministic basis with deterministic coefficients (M 3G)
A Galerkin approach can also be considered for the case that contains undamped deterministic eigenvalues and eigenvectors. For this case, the response vector is defined as follows where γ 0j (ω) and φ 0j correspond to the deterministic scalars and the undamped deterministic eigenvector introduced in Equation (39). c j (ω, θ) ∈ C are unknown constants which need to obtained for each realisation of each frequency. Similarly to the two preceding methods, by applying a sample based Galerkin approach the unknown constants can be computed. Sequentially, the following set of equations is required to be solved for every realisation in each considered frequency where . . n r and c 3 (ω, θ) is the vector that contains the unknown constants c j (ω, θ) The computational effort associated with this method is considerably lower than the other Galerkin methods as the scalars γ 0j only need to be calculated once for each given frequency. The aim of this method is to incorporate the whole stochastic nature of system within the unknown scalars c j (ω, θ). However it is not known if the Galerkin approach can substantially lower the error as all the eigensolutions used are deterministic. This method is of significant interest as it is known that the behaviour of deterministic and stochastic systems can differ substantially especially if the coefficient of variation is significantly large.

Calculation of the response statistics
Calculating response statistics such as the mean and covariance matrices can be deemed funda- The mean of the modulus of the response vector for a given frequency is defined as Similarly the covariance matrix of a given frequency can be defined as follows whereū corresponds to the mean of the response vector and (⋆) T denotes the transpose of ⋆.

Response statistics for the M 1 and M 1G methods
The case of a stochastic projection incorporated with stochastic eigensolutions and a Galerkin error minimisation approach (M 1G) is initially considered. For this case, the response vector can take the following formũ where α j (ω, θ) corresponds to the random scalars seen in Equation (56). Definingũ 1G (ω, θ) in such a way supports neatness whilst calculating response statistics. When analysing the displacement of structures, the modulus of the response corresponds to the maximum displacement or the amplitude seen at each node. The mean vector of the amplitude can be expressed as where |⋆| denotes the absolute value of ⋆. By using the expression from Equation (69), it is possible to define the covariance matrix of the amplitude. In the expression below, rather than using the standard sigma notation, the covariance matrix is defined by C in order to clarify between the covariance matrix and the summations Due to the stochastic nature of this approach, the expression for the covariance matrix cannot be simplified. The mean and covariance matrix of the M 1 method takes a similar form to Equations

Response statistics for the M 2 and M 2G methods
Similarly to the previous case, the method which incorporates the undamped random eigenvalues, undamped deterministic eigenvector and the Galerkin error minimisation approach (M 2G) is initially considered. We will consider the method's response vector in following form where β j (ω, θ) corresponds to the random scalars seen in Equation (62). Due to the undamped deterministic eigenvectors, the mean vector of the amplitude can be expressed as In a similar manner to the previous case, the covariance matrix of the amplitude can be expressed by where the stochastic terms arising in the covariance matrix are given by Due to the undamped eigenvectors being deterministic, less computational effort is needed to compute the expected values arising in the above expression in comparison to that given by Equation (70). The mean and covariance matrix of the M 2 method takes a similar form to Equations (72) and (73). The only difference is that c j = 1 ∀ j = 1, 2, . . . n r and c k = 1 ∀ k = 1, 2, . . . n r .

Response statistics for the M 3 and M 3G methods
When analysing the corresponding mean and covariance matrix arising from the M 3 and M 3G methods, substantial computational reduction can be seen due to the introduction of additional deterministic terms. The response vector of the M 3G method is initially considered in the following where γ 0j (ω) corresponds to the deterministic scalars seen in Equation (64). The mean vector of the amplitude can be expressed as where the only expected value arising is that of the scalars c j (ω, θ). The same is true when considering the covariance matrix of the amplitude where in this instance, the elements of the stochastic term Ψ 3G only contain the scalar terms c j (ω, θ) and c k (ω, θ).
As there are no stochastic terms arising in the response when using the M 3 method, it can be easily deduced thatū 3G (ω) =ũ 3G (ω). Consequently, the covariance matrix would take the form of the zero matrix.  where the values of Z 1 (θ, ω), Z 2 (θ, ω), Z 3 (θ, ω) and y 1 (θ, ω), y 2 (ω), y 3 (ω) are given in Section methods are the Galerkin coefficients c j (θ, ω). It is hoped that by computing these Galerkin coefficients the induced errors will be significantly reduced. Although the values of the Galerkin coefficients need to computed for each θ ∈ Θ and ω ∈ Ω, this size of the linear system which needs to be solved is much lower than the size of the linear system associated with the direct Monte Carlo approach. The computational complexity associated with inversions for θ ∈ Θ and ω ∈ Ω are O(n 3 r ) and O(N 3 ) respectively where n r < N . In the subsequent section, the proposed methods are utilised to analyse two case systems.

Case studies
The six methods proposed in Sections 3 and 6 are applied to two classical structural dynamic systems. The first being an Euler-Bernoulli cantilever beam, and the second being a Kirchhoff-Love plate. Both of the analysed structures have stochastic properties. The stochastic finite element method has been applied in order to discretise both systems. Although two clamped-free systems are examined, the proposed methods can easily be utilised to analyse clamped-clamped systems. This would be executed by adapting the stochastic finite element method.

Euler-Bernoulli cantilever beam
We initially apply the proposed methods to a cantilever beam. Therefore, the displacement and rotational degrees of freedom at the clamped end of the beam are zero. The length of the beam under consideration is 1.00 m and its cross-section is a rectangle of length 0.03 m and height 0.003 m. A harmonic point load is applied at the free tip of the beam. Figure 1 illustrates the configuration. The uncertainty is introduced through the beam's bending rigidity, EI. It is assumed that the bending rigidity is a stationary Gaussian random field of the following form where EI denotes the deterministic value of the bending rigidity. The function a(x, ω) represents a stationary Gaussian field with zero mean, where x is the coordinate direction along the length of the beam. The covariance function is given by where µ a is the correlation length and σ a is the standard deviation. The random field a(x, ω) can be expressed by the Karhunen-Loeve expansion given by Equation (1). The number of terms considered to represent the discretised random field is given by M = 4. Based on the reasoning given in Section 5, for each of the proposed methods only twelve terms have been used in the summations, hence n r = 12. This is a vast reduction as 188 terms have been discarded from each method. For the deterministic case, the distribution of the natural frequencies of the cantilever beam is given in Figure 1. The analytical natural frequencies have been computed by using the methodology given in [57] whilst the finite element solution is computed by solving Equation (42). As to the values obtained by using the finite element method technique match the analytical values, the finite element method technique has been used for the reminder of this study.
The twelve natural frequencies used in the comparison methods are highlighted in Figure 1. As we are only considering the first twelve terms in each of the summations, the methods implementing the sample based Galerkin error minimisation technique requires a linear set of 12 × 12 equations to solved for each sample. The mean vertical amplitude of the displacement at the tip of the beam is shown in Figure 2.
Barring the M 3 method, all the projection methods seem to mimic the results obtained by using the DM CS approach. The M 3 method seems to over estimate the mean vertical amplitude at both the resonance and the antiresonance frequencies of the system. The disparity between the   This frequency value corresponds to the fourth resonance value. All methods seem to mimic that of the DM CS pretty well when σ a = 0.05, however a disparity between the different methods is    The effect of truncating all the projection methods is further explored in Table 2 and Table   3. The approximate L 2 relative error of the mean and standard deviation is explored for different values of n r at a frequency of 42 Hz. This frequency corresponds to the third resonance value. Similarly to the previous case, the bending rigidity has been assumed to be the only stochastic parameter. The bending rigidity of the plate, D, is assumed to be a stationary Gaussian random field of the form D(x, y, θ) = D(1 + a(x, y, θ)) where a(x, y, θ) is a stationary Gaussian field with zero mean and x, y are the coordinate directions of the length and width of the plate. D corresponds to the deterministic value of the bending rigidity of the plate. The correlation function of the random field is assumed to take the following form C a (x 1 , x 2 , y 1 , where σ a is the standard deviation of the bending rigidity, and µ x and µ y are the correlation lengths for both the x and y directions respectively. The forcing vector is again deemed deterministic with a unit norm. This is applied as a harmonic point load at coordinate (0.42, 0.00). The deterministic modal damping matrix consists of a 2% damping factor for each mode. Figure 9 illustrates the configuration of the plate.  The mean and the standard deviation of the vertical deflection amplitude of the plate is further analysed at one of the free corners (0.50, −0.28). This coordinate is labelled P in Figure 9. The mean of the vertical deflection amplitude at point P is illustrated in Figure 10. When σ a is set to 0.05, a good agreement between the DM CS method and all the projection methods is visible.
However when σ a is increased to 0.15, an agreement between the M 3 and DM CS methods is no longer seen. Figure 11 shows the standard deviation of the vertical amplitude of the plate at point P . Due to the M 3 method being deterministic, the standard deviation of this method has not been illustrated.
Barring the M 2 method, a good agreement can be seen between the DM CS method and the other projection methods. Therefore it can be deemed that only the M 1 method or Galerkin methods can incorporate the true standard deviation of the amplitude at point P .
The approximate L 2 relative error of the mean of the response vector is given by Figure 12 for . all the proposed methods. Similarly to the cantilever beam example, the M 3 method is extremely erroneous in comparison to the other methods. This is especially true when σ a = 0.15. The approximate relative error seems to amplify as the frequency increases. If the value of n r were to be increased, it is expected that the relative error at the higher frequencies would decrease.
However increasing the value of n r would increase the computational effort. . for both the relative error of the mean and the standard deviation. Therefore it can be deduced that calculating or approximating the random eigensolutions is not necessary in order to obtain the lowest relative errors for the mean and standard deviation. Combining deterministic eigensolutions and the sample based Galerkin approach is sufficient. . The probability density function of the vertical amplitude of the point P is given in Figure   14. This has been calculated for both values of σ a at 168 Hz. This frequency value corresponds to the 16th deterministic resonance value. Due to the deterministic nature of the M 3 method, its probability density function has been omitted. A decent agreement can be seen between all the methods at the lower σ a value, however a greater discrepancy is seen between the methods at the largest value of σ a . When σ a = 0.15 a good agreement is seen between the Galerkin methods and the DM CS method, however the M 1 and M 2 methods are less successful at mimicking the probability density function of the DM CS method.
Vertical tip deflection (m) Probability density function Probability density function

Summary
A summary of the proposed methods is given below: • A set of hybrid projection methods to calculate the response of stochastic dynamic structural systems has been proposed. The different methods analyse the effect of altering the nature of the coefficients and their associated basis.
• Both the coefficients and their corresponding basis have been computed by utilising the stochastic and deterministic eigensolutions of the structural system.
• The computational effort is reduced by approximating the stochastic eigensolutions and by reducing the modal basis.
• To compensate for the error induced by the proposed hybrid model order reduction technique, a sample based Galerkin error minimisation approach is presented.
• If the sample based Galerkin error minimisation approach is omitted it is necessary for both the coefficients and their associated bases to be stochastic in order to capture an accurate response for a system.
• If the Galerkin error minimisation approach is applied calculating the stochastic eigensolutions is unwarranted. The Galerkin method compensates for the error induced while utilising the baseline eigenmodes. Therefore, a Monte Carlo type sample-based method to evaluate the coefficients and their associate basis can be avoided.
• The application of the Galerkin error minimisation approach in conjunction with projecting onto a deterministic basis with deterministic coefficients (M 3G) produces a level of accuracy comparable to any of the other proposed methods. Our study leads us to suggest that this simple approach has significant potential for analysing stochastic structural systems.

Conclusion
A comprehensive set of hybrid projection methods has been proposed in order to solve a stochastic partial differential equation for structural dynamic systems. Following the implementation of a stochastic finite element method, three projection methods have been developed by utilising the random eigenvalue problem. The first method utilises both random eigenvalues and eigenvectors, the second random eigenvalues and deterministic eigenvectors and the third only uses deterministic eigensolutions. In order to reduce the computational effort associated with each of these methods, the random eigensolutions have been approximated by a first order perturbation and only the dominant projection terms have been retained. Due to the approximations and the reduced modal basis, three additional projection methods have been proposed. These methods utilise a sample based Galerkin error minimisation approach in order to lower the error.
The proposed methods have been applied to a Euler-Bernoulli cantilever beam and to a Kirchhoff-Love plate with stochastic properties described with a random field. It is apparent that if the sample based Galerkin error minimisation approach were not to be implemented the stochastic elastic properties of the random eigenvalues and eigenvectors must be retained in order to capture the stochasticity of the governing equation. The application of the sample based Galerkin method compensates the error incurred when only the baseline eigenmodes are utilised.
Thus the significant computational overhead associated with a Monte Carlo type sample-based basis evaluation is avoided when using the proposed Galerkin error minimisation technique. Future research needs to verify these conclusions for larger real-life structures and different models of uncertainties. One of the main difficulties that needs to be addressed is the elicitation of uncertainties in stochastic parameters for real-life structures. This might require complex hierarchical probabilistic uncertainty descriptors which then must be translated into parametrised coefficient matrices. Additionally, the consideration of non-proportional damping matrices and projection in the space of complex eigenmodes [52] will be a fruitful generalisation of the proposed framework.