## 3 Theoretical Formulation
### Chebyshev Spectral Approximation  
One of the core problems of computational mechanics is the discretization and approximation of structural displacement fields. Another core problem for the beam in Figure 1 is the discretization and approximation of the spatial variation law of the material, cross-section, foundation, and load parameters. A Chebyshev spectral approximation method is used to solve these two problems. This method is described as follows.  
The Chebyshev spectral approximation is based on the Chebyshev polynomials of the first kind. The Chebyshev polynomials are a group of recursive orthogonal polynomials defined by  
\[T_{0}(x)=1,\quad T_{1}(x)=x,\quad T_{k+1}=2 xT_{k}(x)-T_{k-1}(x),\qquad k\geq 1 \tag{9}\]  
where \(k\) is an integer and \(x\in[-1,1]\). For the sake of approximating functions defined on an arbitrary interval \([\ell_{1},\ell_{2}]\), define the scaled Chebyshev polynomials as  
\[\mathcal{T}_{k}(x)=T_{k}(\xi(x)),\quad\xi(x)=\frac{2}{\ell_{2}-\ell_{1}}x- \frac{\ell_{2}+\ell_{1}}{\ell_{2}-\ell_{1}} \tag{10}\]Then, an infinitely differentiable and square-integrable function \(y(x,t)\) in the interval \([\ell_{1},\ell_{2}]\) can be approximated using a truncated Chebyshev series expansion as follows:  
\[y(x,t)=\sum_{k=0}^{N-1}a_{k}(t)\mathcal{T}_{k}(x) \tag{11}\]  
Adopting Gauss-Lobatto sampling,  
\[x_{k}=\frac{\ell_{1}+\ell_{2}}{2}-\frac{\ell_{2}-\ell_{1}}{2}\cos[k\pi/(N-1)], \quad y_{k}=y(x_{k},t) \tag{12}\]  
where \(k=0,1,\ldots,N-1\). There exists a one-to-one mapping between the expansion coefficients vector \(\mathbf{a}=\{a_{k}\}\) and the sampling value vector \(\mathbf{y}=\{y_{k}\}\), which can be expressed as in reference [28]:  
\[\mathbf{a}=\mathbf{\Gamma}_{\text{F}}\mathbf{y},\quad\mathbf{y}=\mathbf{\Gamma}_{\text{B}}\mathbf{a} \tag{13}\]  
where \(\mathbf{\Gamma}_{\text{F}}\) is an \(N\times N\) forward transformation matrix, and \(\mathbf{\Gamma}_{\text{B}}\) is its inverse matrix.  
The \(n\)-th spatial derivative of \(y(x,t)\) can be obtained by the following [28]:  
\[\mathbf{y}^{(n)}=\mathbf{Q}_{n}\mathbf{y} \tag{14}\]  
where \(\mathbf{Q}_{n}\) is the \(n\)-th derivative matrix, a \(N\times N\) matrix with constant elements determined by the interval \([\ell_{1},\ell_{2}]\), and the number of the Chebyshev polynomials used, i.e., \(N\).  
The inner product of two functions \(f(x)\) and \(g(x)\), both of which are infinitely differentiable and square-integrable, can be calculated by the following [29]:  
\[\int_{l_{1}}^{l_{2}}f(x)g(x)\text{d}x=\mathbf{f}^{\text{T}}\mathbf{V}\mathbf{g} \tag{15}\]  
where \(\mathbf{V}\) is called the Chebyshev inner product matrix, a diagonal constant matrix determined by \(N\) and \([\ell_{1},\ell_{2}]\).  
Likewise, the weighted inner product between \(f(x)\) and \(g(x)\) about a weight function \(h(x)\) can be given by the following [29]:  
\[\int_{\ell_{1}}^{\ell_{2}}f(x)h(x)g(x)\text{d}x=\mathbf{f}^{\text{T}}V_{h}\mathbf{g} \tag{16}\]  
where \(\mathbf{V}_{h}\) is the weighted inner product matrix about the weight function \(h(x)\), which is also a diagonal constant matrix determined by \(N\), \([\ell_{1},\ell_{2}]\), and the values of \(h(x)\) at the Gauss-Lobatto sampling points.  
Based on the Chebyshev polynomial approximation and the discrete point sampling in Equations (11) and (12), the calculation formulas of the derivative, inner product, and weighted inner product of the functions in Equations (14)-(16) can be obtained, and the calculus calculation of functions is transformed into matrix operations between discrete vectors. Since \(\mathbf{Q}_{n}\), \(\mathbf{V}\), and \(\mathbf{V}_{h}\) are all constant matrices, the calculation efficiency is greatly improved.