工程上進行CFD計算的時候，常用的離散方法是有限差分法和有限體積法。例如FLUENT采用的就是有限體積法。但是，在需要精度很高的一些基礎研究中，常常采用另一種離散方法——譜方法（spectral method），它具有有限差分法和有限體積法無法比擬的計算精度。今天我們就來聊聊譜方法。

我們以一個具體問題為例。這個問題是常微分方程邊值問題。未知函數u=f(x)滿足

邊界條件是

這個問題是有解析解的，其解析解為

下面我們分別用有限差分法和譜方法來求解這個方程的數值解。讀者可能對有限差分法已經有一定的了解，所以這樣講更有助于讀者對譜方法的理解。

什么是（1）、（2）的數值解呢？（1）、（2）的數值解就是滿足邊界條件（2）的一個函數u*=g(x)，它使得殘差

盡可能地小。不同的數值解法的差別就在于構造數值解u*=g(x)的方法不同，以及評價殘差“小”的方法不同。

首先我們用有限差分法來試試。我們在[-1,1]內劃分n個等距網格節點

其中h=2/(n-1)是網格尺寸。x₁=-1和x_n=1是邊界節點，其余的是內部節點。將各個節點處的數值解u*的值記為u₁,u₂,…,u_n。那相鄰節點之間的區間里面u*的值呢？對于有限差分法來說，其做法是在每個節點附近構造局部的低階插值多項式。例如，在節點x_i附近，讓u*等于一個二階拉格朗日插值多項式p(x)

式子有點復雜，但其實本質很簡單，就是構造一個二次函數p(x)，讓它的函數圖像穿過(x_i-1,u_i-1)、(x_i,u_i)、(x_i+1,u_i+1)三個點。圖1以節點數n=5為例畫出了有限差分法中數值解u*的構造方法示意圖?？梢钥闯霾逯刀囗検降木植啃裕鐇₃附近的u*是用基于x₂、x₃、x₄為插值節點的拉格朗日插值多項式表示的。可以看出u*是分段光滑的函數。

圖1  有限差分法中數值解u*的構造方法示意圖。藍色的圓點代表各個節點處的數值解u*的值u₁,u₂,…,u_n，紅色、綠色和青色的粗線代表局部的低階（這里是二階）插值多項式，而黑色的細線代表數值解u*?？梢钥闯鰑*是分段光滑的函數。

數值解u*構造出來了，那怎樣評價殘差的“小”呢？有限差分法的評價方法比較簡單粗暴，即在每個內部網格節點處令殘差等于零：

下面我們就以采用局部二階插值多項式的有限差分法來具體求解一下（1）、（2）的數值解。將（6）中的p(x)求二階導數，并在x_i處求值，得到數值解u*在節點x_i處的二階導數

因此，殘差在各內部網格節點等于零的條件（公式（7））轉化為下面的差分方程（也叫做差分格式）

此外根據邊界條件（2）有

聯立（9）、（10）可以得到關于u₁,u₂,…,u_n的線性方程組，求解該方程組便可以得到所研究的問題的數值解。由于局部插值多項式是二階的，插值的時候用到三個插值點，所以我們稱之為三點有限差分格式。圖2是節點數n=33時的數值解與精確解的比較。表1是計算域內的最大誤差隨著節點數變化的情況?？梢钥闯觯濣c數大于33之后，每次將網格尺寸h減半都使得誤差減小到原來的1/4左右。理論上，通過泰勒級數展開分析，可以知道差分格式（9）的誤差是O(h²)的，即h→0時，誤差正比于h的平方，現在數值計算結果說明確實是這樣。

圖2  有限差分法節點數n=33時的解答與精確解的比較。采用三點有限差分格式。

 

表1  誤差隨節點數增加的變化；采用三點有限差分格式。

如果我們想加快收斂速度，那么很自然的思路是提高局部插值多項式的階數。例如，如果我們采用局部四階插值多項式，則數值解u*在節點x_i處的二階導數變成

用這個公式重新構造差分格式并算出數值解，計算域內的最大誤差隨著節點數變化的情況如表2所示。由于局部插值多項式是四階的，插值的時候用到五個插值點，所以我們稱之為五點有限差分格式。可以看出，網格節點數大于33之后，每次將網格尺寸h減半都使得誤差減小到原來的1/16左右。理論分析也可知，這種格式的誤差是O(h⁴)的。

 

表2  誤差隨節點數增加的變化；采用五點有限差分格式。

從上面的分析以及數值計算結果都可以看出，局部插值多項式的階數越高，則計算誤差隨網格節點數增加的收斂速度越快。我們自然想到一個極限——能否以計算域的全部網格節點為插值節點構造插值多項式呢？顯然，如果以計算域的全部網格節點為插值節點構造插值多項式，那么數值解u*就不再是分段光滑的函數了，而是全局光滑的函數。這就是譜方法的基本思想——用全局光滑的函數來構造數值解u*。

下面我們就來按照這個思想試一試。我們令數值解u*等于一個函數圖像穿過全部n個點(x₁,u₁),(x₂,u₂),…,(x_n,u_n)的n-1次多項式p(x)，即p(x)滿足

這可以通過拉格朗日插值法構造出來。圖3以節點數節點數n=5為例畫出了這時數值解u*的構造方法示意圖?？梢钥闯霾逯刀囗検绞侨值?，數值解u*是全局光滑的函數。

圖3  譜方法中，數值解u*的構造方法示意圖。藍色的圓點代表各個節點處的數值解u*的值u₁,u₂,…,u_n，綠色粗線代表全局插值多項式，黑色的細線代表數值解u*。

 

這樣，差分格式（9）變成

因為p(x)是全局的插值多項式，所以式中的p’’(x₂)、p’’(x₃)，…，p’’(x_n-1)每個都是全部節點上的數值解u₁,u₂,…,u_n的線性組合，由于表達式比較復雜，這里不一一列出。聯立（13）、（10）可以得到關于u₁,u₂,…,u_n的線性方程組，求解該方程組便可以得到數值解。其計算結果如圖4所示?？梢钥闯?，結果是十分失敗的——隨著節點數n的增加，數值解并沒有越來越接近于精確解，反而是發散了。

（a）n=9

 

（b）n=17

（c）n=33

圖4  采用全局插值多項式的計算結果。

 

這是為什么呢？？其實，數值解u*的構造方法并沒有問題，問題就在于我們對于殘差“小”的定義不再合適了——殘差在一些離散的節點上等于零，并不一定意味著殘差在整個計算域[-1,1]上就一定很小。在有限差分法中，我們使用分段光滑的函數來逼近未知函數，所以殘差也是分段光滑的函數，實踐證明“讓殘差在一些離散的節點上等于零”是可行的。但是現在我們使用全局光滑的函數來逼近未知函數，因此殘差也是全局光滑的函數，實踐證明以前衡量殘差“小”的方法失效了。

那應該怎樣重新定義這個“小”呢？答案就是將殘差R(x)展開為某種正交多項式系的級數形式，然后令展開式中盡可能多的低階項等于零。例如，最常見的就是將R(x)展開為切比雪夫（Chebyshev）級數：

其各項系數的計算公式為

其中T_i(x)是i階切比雪夫多項式。讀者可以類比周期函數傅里葉級數的各項系數的計算公式。

我們知道，對于周期函數來說，光滑函數的傅里葉級數的系數隨著階數的增加迅速衰減。相似地，在[-1,1]上的光滑函數的切比雪夫級數的系數也是隨著階數的增加迅速衰減的（實際上隨著階數的增加是指數衰減的），所以，我們只要強迫展開式（14）中盡可能多的低階項等于零，那么R(x)就可以變得很小很小。注意這里我們充分利用了“殘差是全局光滑的函數”這一性質。如果殘差不是全局光滑的函數（例如在前面的有限差分法中），這一切無從談起。

由于未知量共有n個（u₁,u₂,…,u_n），而邊界條件將其中兩個（u₁,u_n）的值固定了，因此需要n-2個條件才能確定剩余的未知量。因此我們自然地想到，令展開式（14）中的前n-2項的系數等于零：

總結一下，在這種新的對于“小”的定義下，我們的計算方法歸結為：令數值解u*是函數圖像通過全部點(x₁,u₁),(x₂,u₂),…,(x_n,u_n)的n-1次拉格朗日插值多項式，其中u₁和u_n滿足邊界條件（10），剩余的n-2個未知量u₂~u_n-1則通過這樣的方法確定：將u*代入（4）得到殘差的表達式，然后再將殘差的表達式代入方程（16），將積分積出來之后，轉化為關于u₂~u_n-1的線性方程組。這就是譜方法中的Chebyshev-Tau譜方法。當然，Chebyshev-Tau譜方法在具體實施的時候，常常在形式上稍稍改變，即把數值解u*用切比雪夫多項式系展開，將展開式中的系數作為待求解的未知量，而不是用各節點上的數值解u₁,u₂,…,u_n作為待求解的未知量。但這只是形式上的變化，本質上沒有變，所以我們不再贅述。

圖5是用Chebyshev-Tau譜方法，n=17時算出的數值解?？梢钥闯?，數值解與精確解很接近，前面遇到的發散問題得以解決。

圖5  采用Chebyshev-Tau方法的計算結果。n=17。

 

從函數圖像上只能看個大概，下面我們來定量看一下Chebyshev-Tau譜方法的精度。表3給出了計算域內的最大誤差隨著節點數變化的情況。

表3  誤差隨節點數增加的變化；Chebyshev-Tau方法。

可以看出，對于譜方法來說，計算誤差隨網格節點數增加的收斂速度是極快的。當n從33增大到65時，誤差竟然減小到原來的將近1/600。理論上來說，n→∞時，譜方法的計算誤差為O(K^-n)，其中K是一個大于1的常數。這意味著譜方法的誤差是以指數規律減小的，這與有限差分法中誤差以冪函數規律減小的特性是完全不同的，也是任何階的有限差分法都無法比擬的。譜方法的這種收斂特性叫做譜精度（spectral accuracy）。從前面的分析中容易看出，譜方法能達到譜精度的根本原因在于采用了全局光滑的函數來逼近未知函數。采用全局光滑的函數逼近未知函數之后，殘差也是全局光滑的函數，所以殘差的Chebyshev級數的系數隨著階數的增加是指數衰減的。因此，我們強迫級數的前n-2項等于零，其結果就是讓殘差隨著n的增加以指數規律減小。

上面結合具體的問題介紹了Chebyshev-Tau譜方法，這實際上只是譜方法的冰山一角。譜方法的種類其實很多，可以從離散過程以及展開式所用的基函數這兩個方面來分類。從離散過程分類，有Tau方法、Galerkin方法和配置點法三種類型。至于基函數，當計算域是有限區域的時候，經常使用的就是上面介紹的Chebyshev多項式，所形成的譜方法稱為Chebyshev譜方法，但是如果計算域是無限的且具有周期性，則往往選用三角函數來作為基函數，所形成的譜方法稱為傅里葉（Fourier）譜方法。此外還有無限且不具有周期性的區域、半無限區域等等，也有相應的基函數，這里就不討論了。無論哪種譜方法，由于都是用全局光滑的函數來逼近未知函數，所以其計算精度都是達到譜精度的。

譜方法的早期研究有1938年Lanczos的工作以及20世紀60年代Clenshaw, Elliott, Fox等人的工作^[2]。但是，由于譜方法計算量比較大，一直沒有引起人們的興趣。直到1965年Cooley 和Tukey發現了快速傅里葉變換算法，才改變了譜方法的命運。人們發現，對于Fourier譜方法和Chebyshev譜方法，實際計算的時候可以利用快速傅里葉變換算法來大大提高計算效率。于是，20世紀70年代譜方法得以迅速崛起。那時候，美國數學家Steven Alan Orszag是其中最有力的推動者。他在1971年發表于Journal of fluids mechanics的文章“Accurate solution of the Orr-Sommerfeld stability equation”中（文獻[8]），利用Chebyshev譜方法對平面泊肅葉流動的臨界雷諾數進行了精確的計算，通過求解Orr-Sommerfeld特征值問題，得出臨界雷諾數為5772.22的結果。雖然整個研究的結果僅此一個數值而已，但是其影響卻十分深遠。這項工作不僅是譜方法應用于流動穩定性分析的開山之作，同時實際上也確立了譜方法在科學計算中的重要地位。譜方法作為一種求解微分方程的方法，已經廣泛用于流體力學、量子力學、振動、線性和非線性波等等很多領域。下面我們主要談談流體力學方面的應用。

譜方法的特點是精度高，隨著網格加密的收斂速度快；但也有一些明顯的缺點，例如不能適應復雜形狀的計算域等，由于上面討論的例子是一維的，所以看不出這個缺點來。由于這個缺點，譜方法無法應用于實際工程中具有復雜形狀的流體區域。但是在很多流動機理的研究中，計算域簡單，且需要高精度，正適合譜方法施展拳腳。在流體力學中，譜方法最令人矚目的應用是湍流直接數值模擬與流動穩定性分析。

湍流直接數值模擬（Direct numerical simulation，DNS）^[4,5]就是不使用任何湍流模型，而是通過直接求解Navier-Stokes方程組來求解湍流的運動。由于湍流中的旋渦運動的尺度范圍極其寬廣，所以DNS需要很高的網格分辨率，因此DNS自然成為了流體力學中最需要高精度算法的領域之一。湍流直接數值模擬的網格分辨率總是隨著計算機的發展而逐步增加，每個時期從事DNS的研究者都總是希望用上當時最好的計算機。以均勻各向同性湍流為例，最初（1972年）Orszag和Patterson對于均勻各向同性湍流的模擬^[7]是在32³的網格上進行的計算，而21世紀初，對于均勻各向同性湍流的最大規模的計算已經用到了4096³的網格。盡管如此，離模擬工程上的高雷諾數湍流還差很遠。在對物理現象的分辨率不變的前提下，提高算法的精度，就可以采用相對疏一些的網格，從而減小計算時間和內存消耗。

圖6是She等^[6]1990年利用譜方法對均勻各向同性湍流進行直接數值模擬得到的渦結構圖，在該模擬中，計算域為一個正方體，x、y、z每個方向上均采用周期性邊界條件，網格分辨率為96³。由于邊界條件為周期性，所以采用了Fourier譜方法。

圖6  用譜方法進行DNS得到的各向同性湍流中的渦結構（復制自[6]）。該工作發表于1990年Nature。

 

譜方法在流體力學中還常常用于流動穩定性分析。流動穩定性研究的是在Navier-Stokes方程組和某些邊界條件下，數學上適定問題的層流解答，在施加擾動的情形下，擾動衰減還是放大的問題。公眾號之前的文章“層流為何會轉變為湍流：托爾明-施利希廷波的故事”就是這類問題的例子，即平板層流邊界層流動在施加擾動的情況下，擾動最終究竟是會衰減還是會放大。這個問題的重要性在于它揭示了湍流的產生機制，即為什么隨著流動雷諾數的增大，穩定的層流流動會失穩轉化為湍流流動。

流動穩定性分析往往最終歸結為微分方程的特征值問題，例如平行剪切流線性穩定性問題最終歸結為Orr-Sommerfeld特征值問題。對于特征值問題來說，計算精度十分關鍵，這是譜方法在這類問題中得到重要應用的原因。歷史上，Orr-Sommerfeld特征值問題的求解曾經遇到過很大的困難，還吸引過像Heisenberg（海森堡）、Tollmien這樣的大科學家為之努力。其解析求解方法在數學上十分復雜。但是，隨著高速計算機的發展和數值方法的進步，求解Orr-Sommerfeld特征值問題已經不再困難^[9]。今天，譜方法已經稱為求解Orr-Sommerfeld特征值問題的標準方法。

墨爾本大學的研究生劉麗媛和北航宇航學院的研究生胡濤閱讀了本文的初稿并提出了很好的修改意見，在此表示感謝。

 

---

參考文獻

[1]顏慶津. 數值分析. 北京航空航天大學出版社, 2000. 第193至203頁.

[2] Lloyd N. Trefethen. SpectralMethods in Matlab. Oxford University, 2000.

[3] Jie Shen and Tao Tang. Spectral and High-Order Methods with Applications.Science Press, 2006.

[4] M Y Hussaini, T A Zang. SpectralMethods in Fluid Dynamics. Annual Review of Fluid Mechanics, 1987

[5]張兆順, 崔桂香, 許春曉. 湍流理論與模擬. 清華大學出版社, 2005.

[6]Zhen-Su She, Eric Jackson & Steven A. Orszag. Intermittent vortexstructures in homogeneous isotropic turbulence. Nature 344, 226–228 (1990)

[7]Orszag, S. A., Pallerson, G. S. 1972. Numerical simulation of three-dimensionalhomogeneous isotropic turbulence. Phys. Rev. Lett. 28 : 76--79

[8] Orszag,S. A. Accurate solution of the Orr-Sommerfeld stability equation. Journal ofFluid Mechanics, vol. 50, p.689-703

[9]莊禮賢 等. 流體力學. 第二版. 中國科學技術大學出版社, 2009