《數值天氣預報》中非線性計算不穩定的舉例

摘要：數值天氣預報是大氣科學的一個重要分支，是一門實用性很強的應用基礎學科。數值天氣預報本質是用數值方法求解非線性的大氣運動方程組。但當采用格點差分來表示微分方程中的非線性項時，易產生非線性計算不穩定現象。本文以一個簡單的一維非線性平流方程的數值求解過程為例，給出隱式格式差分方程的基本計算方法并重新演示了非線性計算不穩定現象。由于舉例更為簡單，可加深學生對顯式格式和隱式格式優缺點的理解。

關鍵詞：大氣科學；非線性計算不穩定；隱式格式

中圖分類號：G642.0 文獻標志碼：A 文章編號：1674-9324（2013）42-0160-02

一、引言

根據物理定律，比如牛頓第二定律、質量守恒定律、能量守恒定律、氣體試驗定律等，可以得到支配大氣運動的基本方程組。但由于大氣運動基本方程組是一組高度非線性的偏微分方程組，很難求得其解析解，人們可通過數值的方法求其近似解（數值解），這就是數值天氣預報。數值天氣預報是一門實用性很強的應用基礎學科[1]。通過差分方法求解大氣運動基本方程組時，人們發現數值積分過程中會產生計算不穩定問題，這就需要采用恰當的差分格式或積分格式。在《數值天氣預報》課程關于時間積分格式的講解中[2]，重點介紹了“顯示”格式，即差分方程的右端項全部為當前（或和過去）時刻的變量值，通過積數值分可求出方程左端的未來時刻變量值。但對于“隱式”及“半隱式”格式只是簡單提及其概念和特點，比如隱式格式為用未來時刻變量值求出未來時刻變量值，它具有計算穩定、但計算復雜的特點。

在《數值天氣預報》課程中，需講解大量的公式推導和講解，若不能配合簡單而又形象的舉例和圖形，學生尤其是本科生作為授課對象，將很難理解和接受本課程中的相關內容，講課的效果也將大打折扣。在非線性不穩定計算的舉例中，教材中[2]雖給出了采用不同的初值和不同的差分方案對計算穩定性的影響，但并沒有清楚地列出其求解過程，因此學生很難了解隱式格式差分的具體求解過程，對教材中列出的“顯示”和“隱式”格式的各自優缺點更是難以理解。因此，需要對兩種格式的計算過程進行相應的講解，尤其是隱式格式。此外，教材中[2]對同一個微分方程構造的兩個不同的差分方程中，除隱式格式和顯式格式的差異外，還存在著對■采用了不同的差分格式，即顯式格式采用中央差格式，隱式格式采用前差格式。本課程[2]已清楚的講解到中央差格式雖具有較高的計算精度，但在時間差分計算時存在計算解的問題，若初值取得不當，則計算解會有較大的振幅。因此，從邏輯上講，教材中給出的不同的差分方案的影響，實際上不僅僅來源于顯式格式和隱式格式的差異，還來源于對時間微分采用不同差分格式的差異。這又加大了學生對顯式格式和隱式格式特點的理解難度。

針對上述問題，本文將以簡單的一維非線性平流方程為例，給出隱式格式差分方程的具體求解過程，重新探討非線性計算不穩定現象，目的是使學生更好地了解顯式格式和隱式格式差分方程的求解過程，深刻理解兩種格式各自的優缺點。

二、非線性計算不穩定的計算實例

以大氣科學中極具代表性的一維平流方程為例：■+u■=0，0≤x<1。它可以寫成

■+■（■）=0，0≤x<1 （1）

或■+■（u■+■），0≤x<1 （2）

以上兩式與教材[2]基本一致，不同的是這里的x取值范圍并不到1。在大氣科學中，方程或模式的計算可在全球或某一緯圈上進行。在該情況下，沒有緯向側邊界條件。對于上式而言，可認為u在x=1的取值等于u在x=0的取值，也即循環邊界條件。在構造上述微分方程相應的差分方程過程中，對（1）式和（2）式分別采用顯式格式和隱式格式：

uin+1=uin-■[（ui+1n+uin）2-（uin+ui-+1n）2] （3）

uin+1=uin-■[（■i+1+■i+■i-1）（■i+1-■i-1）] （4）

其中上標n為第n步，下標i為第i個格點，■i=（uin+1+uin）/2。可見，與教材中不同的是，（1）式和（2）式中■均取了前差格式，這樣可避免由于三個時間層計算而出現的計算解問題，有利于問題的討論更加集中。

同樣給定兩種不同的初值，兩者僅相差一個常數：

ui0=sin2πiΔx （5）

ui0=1.5+sin2πiΔx （6）

計算中，Δx取=1/3，Δt=0.004，則|u■|≤umax=|u■|=2.5×0.004×3=0.003<1，根據線性穩定性判據，其滿足線性穩定性條件。將兩種差分方案與兩種初值兩兩組合，可得四種計算結果。在線性穩定性條件情況下，若結果仍出現不穩定現象，則一定是非線性計算不穩定。為檢查計算穩定性，計算每一步中所有格點的總動能■∑ui2，若該值在某一有限區域內變化，則計算穩定，否則不穩定。

三、隱式格式的求解

顯式差分方程（3）的求解過程即是將已知的n時刻u值代入等式右端算出等式左端未知的n+1時刻u值，可見，求解過程簡單。至于隱式差分方程（4），其求解過程，較復雜。首先將（4）式寫在[0，1）的x0=0、x1=1/3和x2=2/3三個格點上，并令m=-■，u0n+1+u0n=X，u1n+1+u1n=Y，u2n+1+u2n=Z，u0n=a，u1n=b，u2n=c。可見，a、b、c均為已知的第n步值。采用循環邊界條件可得三元二次方程組：

X=2a+m（X+Y+Z）（Y-Z）Y=2b+m（X+Y+Z）（Z-X）Z=2c+m（X+Y+Z）（X-Y） （7）

將（7）式中的三式相加可得：X+Y+Z=2a+2b+2c，再令m（2a+2b+2c）=d，該d值也是已知的第n步值，（7）式可化為三元一次方程組：

X=2a+d（Y-Z）Y=2b+d（Z-X）Z=2c+d（X-Y） （8）

最終可利用已知的a、b、c和d值分別求得n+1步未知的Z、Y、X值：

Z=■Y=■X=2a+d（Y-Z）（9）

再分別將其減去c、b和a值，可得n+1步的u2n+1、u1n+1和u0n+1。由此可見，隱式差分方程的求解過程較為復雜。需指出的是，本文在[0，1）僅選取了3個格點，若選取教材中的10個點（Δx=0.1），則需在10個格點上寫出10個差分方程，并進行聯立，求解十元一次方程組，其求解過程更為復雜。

四、計算結果及分析

圖1a和1b分別給出了兩個初值、兩種計算方案的計算結果。初值取（5）式用顯式方案（3）式的計算結果表明（圖1a實線），動能逐步增大，在500步以后突然急劇增加，出現按指數增加的趨勢；但若給初值加上一個常數后（圖1b中實線），總動能在4m2/s2左右變化，表明計算結果穩定。至于隱式方案（4）式，無論取哪種初值，結果均穩定。

至于產生如圖1a中的不穩定現象，仍可利用混淆誤差理論進行解釋，即網格系不能正確分辨短波長的波動而導致不穩定。本文例子在[0，1）的一個周期范圍內僅取三個格點，采用循環邊界條件，即將[0，1）進行I=3等分。可見，該網格系只能正確識別平均值0波、波長為的3/2波和波長2Δx為3Δx的1波波動。若波數k1=k2=3/2的兩個波動相互作用，則可產生0波和3波的波動。其中3波波動超出該網格系的識別能力，將會被錯誤的識別為0波。該0波的能量將不斷的積累，從而可導致不穩定現象。該過程也可通過（10）式得到驗證：

sin■=sin■=0=sin■cos■=cos■=1=sin■ （10）

因此，雖然本文的計算結果與教材中基本一致，但舉例十分簡單，這有利于學生的理解和接受。

五、結束語

本文通過簡單的計算實例重新探討了差分格式對非線性計算穩定性的影響。這里的“簡單”，主要指將一維非線性平流方程的時間偏導項統一地取成前差格式，同時，差分方程僅寫在三個格點上。從而，隱式格式差分方程的求解過程便成為三元一次方程組的求解過程。該求解過程比顯式格式差分方程復雜，但計算結果穩定，充分體現出隱式格式和顯式格式的優缺點。雖最終的計算結果與教材中[2]基本一致，但本文舉例更為簡單、易懂，且給出了詳細的求解過程，有助于學生自己動手推導和計算求解，以加深其對顯式和隱式格式的理解，并深刻體會各自的優缺點。

參考文獻：

[1]紀立人.數值天氣預報發展進程中若干亮點的回顧及其啟迪[J].氣象科技進展.2011，1（1）：40-42.

[2]沈桐立，田永祥，葛孝貞，陸維松，陳德輝.數值天氣預報[M].北京：氣象出版社，2003.9.