# Structural Reliability Analysis Using Stochastic Finite Element Method Based on Krylov Subspace
## 2 Theoretical Development  
For a structure with \(n\) degrees of freedom (DOFs), its displacement under static load is usually obtained by solving the following linear equations:  
\[K\cdot x=l \tag{1}\]  
where \(K\) denotes the original stiffness matrix of the structure, \(x\) denotes the original displacement vector, and \(l\) denotes the static load vector. From Equation (1), the original displacement vector can be calculated using  
\[x=K^{-1}l \tag{2}\]  
As mentioned earlier, many physical quantities in the structural finite element model, such as the elastic modulus of materials, are random variables. This leads to the stiffness matrix of the structure also being a random matrix. The structural stiffness matrix considering randomness can be expressed as the sum of the stiffness matrix of each element multiplied by a random coefficient, that is  
\[K_{d}=\sum_{i=1}^{N}c_{i}K_{i} \tag{3}\]  
where \(K_{d}\) denotes the random stiffness matrix, and \(c_{i}\) represents a random coefficient that obeys a certain distribution, such as a normal distribution and a uniform distribution. For a random FEM, its displacement under a static load is also obtained by solving the following linear equations as  
\[K_{d}\cdot x_{d}=l \tag{4}\]  
where \(x_{d}\) represents the displacement of a random FEM sample under the static load. Equation (4) can still be solved using the conventional algorithm as:  
\[x_{d}=K_{d}^{-1}l \tag{5}\]  
The approach of calculating the structural displacement under each random sampling through the direct inversion of Equation (5) is called the complete analysis method, which is also known as the classical Monte Carlo finite element method. The result obtained using this method is the exact solution, which is usually used to validate the computational accuracy of other stochastic analysis methods. However, for large-scale engineering structures, if the displacement of each random FEM sample is calculated using Equation (5), the calculation cost will be very expensive. In order to save on computation, some fast stochastic finite element methods have been developed in the past decades, such as the Neumann stochastic finite element algorithm (NSFE). The basic idea of the NSFE method is to use Neumann series expansion to quickly obtain the approximate solution of \(K_{d}^{-1}\). To this end, \(K_{d}\) is rewritten as  
\[K_{d}=K+\Delta K \tag{6}\]where \(\Delta K\) denotes the change in the stiffness matrix. Substituting Equation (6) into (5) yields  
\[x_{d}=(K+\Delta K)^{-1}l \tag{7}\]  
Equation (7) can be expanded by using the Neumann series [27; 28] as  
\[x_{d}=[K^{-1}-(K^{-1}\Delta K)K^{-1}+(K^{-1}\Delta K)^{2}K^{-1}-\cdots]\cdot l \tag{8}\]  
Using Equation (8), the NSFE solutions with different accuracies can be obtained by retaining the first few terms. For example, the NSFE solution with \(m\) items can be expressed as:  
\[x_{d}\approx[K^{-1}-(K^{-1}\Delta K)K^{-1}+\cdots+(K^{-1}\Delta K)^{m-1}K^{-1} ]\cdot l \tag{9}\]  
One can conclude from Equation (9) that the NSFE method avoids directly inverting the random stiffness matrix \(K_{d}\) through the Neumann series, thus saving the calculation amount. As noted in the literature [22; 29; 30; 31; 32], the computational efficiency and accuracy of the NSFE method are primarily contingent upon the spectral radius of the perturbation matrix \(\Delta K\). Specifically, when the spectral radius of \(\Delta K\) is small, a relatively modest number of terms \(m\) in the NSFE calculations suffices to yield highly precise results. Conversely, if the spectral radius of \(\Delta K\) is substantial, a significantly larger number of terms \(m\) is required in the NSFE computations to achieve comparable levels of precision. Consequently, the NSFE method is frequently regarded as unsuitable for tackling problems characterized by significant variance in random variables, referred to as \"large random fluctuations\".