一、“近似”的兩種分類
一個復雜的函數，可以通過一系列的“基底函數”（base function）的組合來近似，也就是函數逼近，有兩種典型的方法：
- 基於全域的逼近，如傅立葉級數展開；
- 基於子域的分段函數組合，如有限元方法。
第一種函數逼近方式，就是力學分析中經典的瑞利-里茲方法（Rayleigh-Ritz），這種方法的特點是基底函數比較復雜，一般是高階連續函數，通常僅需采用前面幾階函數組合即可得到較高的逼近精度，比如展開為傅立葉級數。
第二種函數逼近方式就是現代力學分析中的有限元思想，即“分段逼近”，每一個分段函數一般比較簡單，采用線性函數或者二次函數即可，但是需要較多的分段才能得到逼近效果，工作量比較大。
關於兩種逼近方法更加形象的說明，可以參考曾攀老師書中的一個圖片，也就是下面這張圖
1.1 利用基於全域的逼近方法求近似解
瑞利-里茲法的核心觀點就是選取“試函數”，這個試函數必須首先滿足位移邊界條件，當然還會帶有一些待定系數，然后將其他的變量都用這個“試函數”來表達，通過其他的邊界條件或者能量方法（虛功原理或極小勢能原理）來求解待定系數。
1.1.1 利用虛功原理求解近似解
簡單來說，虛功原理的含義就是如果系統有一個虛位移，那么外力在虛位移上所做的外虛功應該內力所做的內虛功。
內虛功
\begin{equation}
\delta U=\intop_{\varOmega}\sigma_{x}\delta\varepsilon_{x}d\varOmega
\end{equation}
外力虛功
\begin{equation}
\delta W=\intop_{l}p\delta vdx
\end{equation}
由於選取的“試函數”只需滿足位移邊界條件即可，所以“試函數”的選取條件非常寬泛。但是不同的試函數求得最終結果的精度確實相差很大的，所以選取一個合理的試函數非常關鍵。
比如，對於純彎曲的簡支梁而言，它的位移邊界條件就是在0和L處的位移為零，而且位移的形狀應該是中間大，然后逐漸向兩端減小為0。很自然的首先想到的就是二次函數拋物線，另外一種就是三角函數正弦曲線。為簡化，這里將梁的長度L均取為1
- 試函數$v1(x)=-x^{2}+x$
- 試函數$v2(x)=\frac{1}{4}\sin(\pi x)$
畫出上述兩個試函數在0~1之間的圖形
可以看到二者都符合位移邊界條件，而且形狀的大概模樣和預想的梁撓度曲線是一致的。
當試函數采用$c_{1}\sin\left(\frac{\pi x}{L}\right)$時，根據虛功原理，可以求得系數c1，得到最終的撓曲線為
\begin{equation}
\frac{4L^{4}p_{0}}{\pi^{5}\text{EI}}\sin\left(\frac{\pi x}{L}\right)
\end{equation}
當試函數采用$c_1 \left(L x-x^2\right)$時，得到的最終撓曲線為
\begin{equation}\frac{L^{2}p_{0}}{24\text{EI}}(Lx-x^{2})
\end{equation}
令$\frac{p_{0}}{EI}=1$和$L=1$，分別畫出0~L直接解析解、采用兩種試函數得到的近似解的圖形
可以看出采用$c_{1}\sin\left(\frac{\pi x}{L}\right)$作為試函數，最終的到曲線v1(x)與解析解v0(x)是非常接近的，在上圖中幾乎看不出差別，但是采用$c_1 \left(L x-x^2\right)$作為試函數求解出的近似解卻和解析解相差比較大，進一步放大圖形
可以看到，當圖形放大之后，才能看出v0(x)和v1(x)直接的細微差別。
1.1.2 利用極小勢能原理求近似解
通過上面采用虛功原理求近似解，可以看出試函數采用三角函數得到的結果更加准確，事實上這就是是傅立葉級數的展開式。
極小勢能原理，也就是勢能駐值原理，簡而言之就是系統的勢能極小時，系統是穩定的。
我們一般求系統的勢能是用應變能減去外力功，即$\varPi=U-W$。
其中U表示應變能，它的表達式為
\begin{equation}
U=\dfrac{1}{2}\intop_{\varOmega}\sigma_{x}\varepsilon_{x}d\varOmega
\end{equation}
W表示外力功，表達式為
\begin{equation}
W=\intop_{l}pvdx
\end{equation}
基於上面虛功原理，可以看出試函數選取三角函數更加合理，事實上，由於解析解是采用冪級數形式的方程，采用傅立葉級數（三角函數）進行近似更加合理，后面我們也會看到，實際上試函數的待定系數就是解析解的傅立葉系數。
選取試函數為
\begin{equation}
v(x)=c_{1}sin\dfrac{\pi x}{l}+c_{2}sin\dfrac{2\pi x}{l}+c_{3}sin\dfrac{3\pi x}{l}
\end{equation}
采用Mathematica進行求解過程如下
\[
\text{v3}(\text{x$\_$})\text{:=}c_1 \sin \left(\frac{\pi x}{L}\right)+c_2 \sin \left(\frac{2 \pi x}{L}\right)+c_3 \sin \left(\frac{3 \pi x}{L}\right)
\]
\[
\text{U3}=\frac{1}{2} \text{EI} \int_0^L \text{v3}''(x)^2 \, dx=\frac{\pi ^4 \left(c_1^2+16 c_2^2+81 c_3^2\right) \text{EI}}{4 L^3}
\]
\[
\text{W3}=p_0 \int_0^L \text{v3}(x) \, dx=\frac{2 \left(3 c_1+c_3\right) L p_0}{3 \pi }
\]
\begin{equation}
\Pi =\text{U3}-\text{W3}=\frac{\pi ^4 \left(c_1^2+16 c_2^2+81 c_3^2\right)
\text{EI}}{4 L^3}-\frac{2 \left(3 c_1+c_3\right) L p_0}{3 \pi }
\end{equation}
將勢能表達式對c1,c2,c3求駐值，即對各自的偏導為0，可以求出c1，c2，c3。
\[
\text{Solve}\left[\left\{\frac{\partial \Pi }{\partial c_1}=0,\frac{\partial \Pi }{\partial c_2}=0,\frac{\partial \Pi }{\partial c_3}=0\right\},\left\{c_1,c_2,c_3\right\}\right]
\]
\[
\left\{\left\{c_1\to \frac{4 L^4 p_0}{\pi ^5 \text{EI}},c_2\to 0,c_3\to \frac{4 L^4 p_0}{243 \pi ^5 \text{EI}}\right\}\right\}
\]
可以看到c1的值與使用虛功原理求得的值是一樣的。由於采用了兩項三角函數組成的試函數，最終求得的撓曲線相比虛功原理一項三角函數得到的撓曲線更加精確。根本的原因是這樣的，上述c1，c2，c3就是解析解的按照傅立葉級數展開之后的傅立葉系數。
根據周期為$2l$的周期函數的傅立葉級數公式（《高等數學》同濟六版下冊P316）
由於解析解為
\begin{equation}v(x)=\dfrac{p_{0}}{EI}(\dfrac{1}{24}x^{4}-\dfrac{1}{12}lx^{3}+\dfrac{1}{24}l^{3}x)
\end{equation}
為了簡便，將$\dfrac{p_{0}}{EI}$都看作1，得到簡化后的撓曲線解析解為
\begin{equation}
v(x)=\dfrac{1}{24}x^{4}-\dfrac{1}{12}lx^{3}+\dfrac{1}{24}l^{3}x\end{equation}
這個函數在[0,1]的圖形可以做奇函數的周期拓展，所以傅立葉級數展開式只有正弦項，沒有余弦項。其傅立葉系數為
\[
v(\text{x$\_$})\text{:=}\frac{L^3 x}{24}-\frac{L x^3}{12}+\frac{x^4}{24}
\]
\[
\text{b1}=\frac{2 \int_0^L v(x) \sin \left(\frac{\pi x}{L}\right) \, dx}{L}=\frac{4 L^4}{\pi ^5}
\]
\[
\text{b2}=\frac{2 \int_0^L v(x) \sin \left(\frac{2 \pi x}{L}\right) \, dx}{L}=0
\]
\[
\text{b3}=\frac{2 \int_0^L v(x) \sin \left(\frac{3 \pi x}{L}\right) \, dx}{L}=\frac{4 L^4}{243 \pi ^5}
\]
可以看出根據傅立葉公式求出的傅立葉系數與根據勢能原理求出的試函數的待定系數是一樣一樣的。
下面一篇會從有限元方法推導梁單元。To be continue………….