## 3 Generation of Cross-Correlated Random Fields
### Methodology
#### 3.1.2 Generation of Random Fields by Numerical Methods  
Initially, random fields are continuous functions with respect to the spatial coordinates. However, to be able to utilize this concept within the framework of material modeling and numerical simulations by the finite element method (FEM), a discrete representation is necessary. Over the years, many different methods have been developed for this purpose. Examples are the midpoint method [37; 38], spatial averaging [39], and the shape function method [40; 41; 42]. Commonly used approaches are based on the series expansion technique [43]. [PERSON] et al. [44] extended the approach by expressing random fields by a spatial decomposition of the correlation functions derived by the well-known Karhunen-Loeve expansion (KLE) [45]. The spatial decomposition obtained by the KLE for a mean free random field with a standard deviation equal to one was given by [PERSON] et al. [46]:  
\[Z(\omega,\mathbf{x})=\sum_{n=0}^{\infty}\sqrt{\lambda_{n}}\phi_{n}(\mathbf{x} )Z_{n}(\omega). \tag{2}\]  
Here, \(Z_{n}(\omega)\) are uncorrelated, standard normally distributed random variables and \(\lambda_{n}\) and \(\phi_{n}(\mathbf{x})\) are the eigenvalues and eigenfunctions of the correlation function kernel, which can be obtained by solving a Fredholm integral equation of the second kind. Since the integral equation:  
\[\int_{\Omega}\mathrm{Cov}[\mathbf{x}^{\prime},\mathbf{x}]\phi_{n}(\mathbf{x}) \,\mathrm{d}\,\mathbf{x}=\sum_{n=1}^{\infty}\lambda_{n}\phi_{n}(\mathbf{x}^{ \prime}) \tag{3}\]  
has closed solutions only for a few types of correlation functions defined on a rectangular domain [47], e.g., the exponential correlation function [31; 48], numerical methods are often required. This is especially the case for multidimensional fields [49]. By deploying numerical integration methods, the eigenfunctions \(\phi_{n}(\mathbf{x})\) as a solution of the Fredholm integral are approximated by a set of functions \(h_{i}\):  
\[\phi_{n}(\mathbf{x}^{\prime})\approx\hat{\phi}_{n}(\mathbf{x}^{\prime})=\sum_ {i=1}^{N}d_{i}^{n}h_{i}(\mathbf{x}^{\prime}), \tag{4}\]  
where the parameters \(d_{i}^{n}\) are unknown and need to be determined. One technique here is the expansion optimal linear estimate (EOLE) developed by [PERSON] et al. [49]. It is based on the linear estimation theory and belongs to the group of series expansion methods. Furthermore, it can be shown that the EOLE is equivalent to the Nystrom method [50] with uniformly distributed integration points [47]. Within the Nystrom method, the integral eigenvalue problem is written as  
\[\sum_{i=1}^{N}w_{i}\operatorname{Cov}[\mathbf{x}_{n},\mathbf{x}_{i}]\hat{\phi}_ {j}(\mathbf{x}_{i})=\hat{\lambda}_{j}\hat{\phi}_{j}(\mathbf{x}_{n}), \tag{5}\]  
which can be reorganized in matrix notation as  
\[\mathbf{C}\mathbf{W}\mathbf{y}_{j}=\hat{\lambda}_{j}\mathbf{y}_{j}. \tag{6}\]  
The \(N\times N\) matrix \(\mathbf{C}\) is symmetric positive semi-definite and holds the elements \(c_{n,i}=\operatorname{Cov}[\mathbf{x}_{n},\mathbf{x}_{i}]\); the integration weights \(w_{i}\) are stored in the diagonal matrix \(\mathbf{W}\) of the size \(N\times N\), and \(\mathbf{y}_{j}\) is a column matrix with the entries \(y_{j,n}=\hat{\phi}_{j}(\mathbf{x}_{n})\). Assuming a uniform distribution of the points \(\mathbf{x}_{i}\) over the domain \(\Omega\) and a equispaced structured grid, respectively, all integration weights \(w_{i}\) are the same. In this case, the matrix \(\mathbf{W}\) reads  
\[\mathbf{W}=w\mathbf{I}, \tag{7}\]  
with \(\mathbf{I}\) being the identity matrix and \(w=|\Omega|/N\). This leads to  
\[\mathbf{C}\mathbf{y}_{j}=\hat{\lambda}_{j}^{*}\mathbf{y}_{j}, \tag{8}\]  
giving the equivalent of Equation (6) for the EOLE, with the eigenvalues \(\hat{\lambda}_{j}^{*}\) and eigenfunctions \(\mathbf{y}_{j}\) of the covariance matrix \(\mathbf{C}\).