1、CH3 數值計算方案（1）譜方法；（2）有限差分方法；1 差分方法和差分格式2 差分格式的基本性質3 時間積分方法和積分格式4 有限差分格式的誤差分析5 非線性計算穩定性主要內容1 差分方法和差分格式首先要對解域進行離散化（包括空間和時間的離散化），建立對應的網格系。一、最簡單的一階線性偏微分方程的數值解法一維線性平流方程：(3.1)0 xuctu 在x-t平面上建立以 x和 t為間隔的網格系，則任意一點的坐標為， 則tnxiutxuunini,NntntMixixni1 , 01 , 0 xtixnt x 為空間格距（步長）， t 為時間格距（步長）。in下面用Taylor展開式來構造出常用

2、的幾種差分格式，據Tailor展開式，當 x為一小量時，有以下關系式：（3.3a）（3.3b）.! 3)(! 2)(),(),(333222xxuxxuxxutxutxxu.! 3)(! 2)(),(),(333222xxuxxuxxutxutxxu二、差分格式的建立（3.3a）Rxx, t ux, tx uxutx,第一種差分格式向前差格式（前差格式）xx+xx-x.! 3)(! 2)(),(),(333222xxuxxuxxutxutxxuxoxxuxxuR.! 3)(! 223322R為截斷誤差，它在一定程度上代表了差分格式的精度。（3.3b）Rxx, tx ux, t uxutx,第二

3、種差分格式向后差格式（后差格式）xx+xx-xxoxxuxxuR.! 3)(! 223322.! 3)(! 2)(),(),(333222xxuxxuxxutxutxxu（3.3a）（3.3b）Rxx, tx ux, tx uxutx2,第三種差分格式中央差格式xx+xx-x.! 3)(! 2)(),(),(333222xxuxxuxxutxutxxu.! 3)(! 2)(),(),(333222xxuxxuxxutxutxxu2233.! 3)(xoxxuR前差格式： R=O x 3.3c 后差格式： R=O x 3.3d 中央差格式： R=Ox2 3.3e R的階次越高，則誤差越小，差分格

4、式的精度越高。中央差格式的精度要比前差格式和后差格式高。一階微分的三種常用的差分格式：Rxx, t ux, tx uxutx,Rxx, tx ux, t uxutx,Rxx, tx ux, tx uxutx2, 同樣，可以用時間差分格式來表示時間微商，得到類似的時間差分格式。這樣，若時間t取前差格式，空間x分別取前差、后差和中央差格式，則可以構造出以下三種差分方程：011xuuctuunininini O t , x 3.13O t , x 3.14O t , x23.15 011xuuctuunininini02111xuuctuunininininiutxu),(niutxxu1),(ni

5、utxxu1),(1),(niuttxu1),(niuttxu二階微分的差分問題：.! 3)(! 2)(),(),(333222xxuxxuxxutxutxxu.! 3)(! 2)(),(),(333222xxuxxuxxutxutxxuRxuuuxuRxx, tx ux, tux, tx uxuninininitx211222,222)(222222yAxAAxy(i,j)(i-1,j)(i+1,j)(i,j+1)(i,j-1)2, 1, 1,22)(2s AAAxAjijijiji21,1,22)(2s AAAyAjijijiji2, 1, 11,1,2)(4)(sAAAAAAjijiji

6、jijiji差分格式一旦建立，這種格式是否是可用的？必須在采用差分格式之前對之進行進一步的考查，討論差分格式的相容性收斂性穩定性.等問題。2 差分格式的基本性質一、差分格式的基本性質當空間步長 x 和時間步長t 很小時，差分方程是否逼近微分方程，這就是差分格式的相容性（一致性）問題。在實際問題中，通常利用R來考察差分格式的相容性（一致性）。如果當t , x 0時，R 0的，則差分方程與微分方程是相容的或是一致的。1、相容性（一致性）問題可見，以上三種差分格式構造的差分方程與微分方程（3.1）是相容（一致的）。在設計一種差分格式時，首先必須考慮其相容性（一致性）。011xuuctuuninini

7、ni O t , x 3.13O t , x 3.14O t , x23.15 011xuuctuunininini02111xuuctuunininini在一定的定解條件下，差分方程的解 是否逼近微分方程的解 f 的問題，稱之為差分格式的收斂性問題。即 必須注意，差分格式的相容性并不能保證其收斂性。由于對應微分方程的真解通常是并不知道的，要直接證明差分格式的收斂性是存在一定困難的，在考慮說明收斂性問題之前，先來討論差分格式的穩定性問題。2、收斂性問題設f 代表微分方程的解， 代表差分方程的準確解， 代表差分方程的近似解（數值解）。那么：解的截斷誤差 舍入誤差FF)()(FFFfFffFxt0

8、,limF在時間積分過程中，由于舍入誤差的影響，差分解的誤差是否隨時間增長的問題，即差分格式的計算穩定性問題。也就是說，當時間步長趨向于0時，在整個求解區域內，舍入誤差是否保持有界的問題，若保持有界則是穩定的。3、穩定性問題舍入誤差有界差分格式的相容性、收斂性及穩定性的關系如何？下面將給出拉克斯（Lax）等價定理：如果差分方程逼近微分方程，即差分格式與微分方程是相容的，或者差分格式滿足相容性條件，差分格式的穩定性，保證了其收斂性（計算穩定性是收斂性的充分必要條件）。 可見，在差分格式滿足相容條件時，其穩定性就保證了其收斂性，因此，要保證差分格式的穩定性也就顯的格外重要了。關于差分格式的穩定性的

9、討論我們將在稍后的內容進行介紹。另外，在設計差分格式時，還必須考慮計算的有效性和節約計算資源，應盡可能選擇計算簡單、速度快、使用計算機內存少的方案。 IkxAexFxu0 ,ctxIkAectxFtxu,0 xuctu 二、差分格式的線性計算穩定性問題一維線性平流方程設其具有波動形式的初值則微分方程的解（3.32）時間用前差分，空間用向后差分格式，可得差分方程O t , x 3.14011xuuctuunininini令于是有：xtc nininiuuu111 ninininiuxtcuxtcuu11而差分方程的近似解，可以表示為：iIkxnnieAu 011xuuctuunininini I

10、kxctxIketAAectxFtxu,把近似解代入差分方程，進一步有：故最終可以將上式寫成：（3.34）令上式可寫成：（3.35）xIknnnIkxnIkxnIkxneAAAeAeAeAiii11111nxIknAeA11GexIk 1nnGAA1根據以上遞推式子，差分近似解可表示為： （3.37）可見，如果上式滿足 ，即 時，則有差分近似解 的振幅將不隨n的增加而增大，即差分近似解是有界的，因此此時的差分格式就是穩定的。這就是馮-紐曼的穩定性必要條件。定義G為增幅因子：不難得出：（3.36）xIknneAAG 1101AGGAAnnniIkxnnieGAu01G11nnAAniu現根據以上

11、討論的知識，用馮-紐曼方法來證明差分格式（3.14）的穩定性。 歐拉公式： 首先，對于該差分格式，其增幅因子為：xIknneAAG 11xixexixeixixsincossincos利用歐拉公式將上式展開：222sincos1sincos1xkxkGxkIxkG 2sin141cos1121cos12121cos12221cos1221sincos12cos1222222222xkxkxkxkxkxkxkxk12G01014 如果要 ，則應滿足： 對上述不等式進行討論：（1） （2）可見，當 時，有 成立，則在這個條件下，差分格式（3.14）是穩定的。 10010 010 101G綜上所述，

12、用Von-Neumann穩定性判別方法來證明差分格式的線性計算穩定性時的主要步驟為：1設解的波動形式，代入差分方程。2得出其對應的增幅因子G。3討論 時的情況。4根據滿足 的條件判斷滿足差分格式穩定性的條件。 1G1G3 時間積分方法和積分格式本節主要介紹時間積分格式的一般性概念，并給出不同的時間積分格式及其穩定性的分析，以便對數值預報模式設計的算法有一個初步的了解。 一、時間積分方法Fdtdu數值預報的模式方程，通常可寫成：（3.55a）此方程為振動類方程，其中 ，c為相速。若令 ,則平流方程就轉化為方程（3.55a）的形式了。0 xActA IkxetuAuIIkcudtdu kc 對于前

13、面討論的一維平流方程：（3.55c）假定可將平流方程化為如下形式：（3.55d）uIF 以振動方程 為例來說明問題。假設 為已知，來構造時間積分格式，來求下一時刻的值 。其中 表示u在 時刻的值， 類推。二、時間積分格式Fdtdu1,nnuu1nunutnt11,nnuuuIF 1、二時間層的積分格式（非迭代格式）2、二時間層的積分格式（迭代格式）3、三個時間層的積分格式主要介紹三類時間積分格式，并討論每種格式的性質：此類時間積分格式僅涉及到兩個時間層，故稱為二時間層的積分格式。1、二時間層的積分格式（非迭代格式）nnntFuu111nnntFuua、歐拉格式：b、后差格式（隱式格式）（3.6

14、3）（3.64）c、梯形格式（隱式格式）（3.65）112nnnnFFtuu其中A+B=1，對于 （a）有：A=1，B=0（b）有：A=0，B=1 （c）有：A=1/2，B=1/211nnnnBFAFtuu上述三種格式可寫成以下通用形式：（3.66）格式分析：即利用Von-Neumann方法，討論其穩定性。（其中 ，進而有 ）uIF nnuIF 11nnnnBFAFtuu根據G的表達式，對時間積分格式的穩定性進行分析。 求其增幅因子為：(3.69）22211111tBtItABtIBtIAuuGnn （3.66）a、歐拉格式：A=1，B=0，故其為絕對不穩定格式。 b、后差格式（隱式格式）：A

15、=0，B=1可見其為無條件穩定格式，或絕對穩定格式。 c、梯形格式（隱式格式）：A=1/2，B=1/2 故為中性格式，其數值解振幅不變。11122tGtIG 11111222tGttIG 1411411222GttItG 2、二時間層的積分格式（迭代格式）1*,11*,nnnnnntFuutFuunnFu nF1*, nu1*, nF1nua、歐拉后差格式（3.72）是一個不斷迭代的過程。具體說明：1*,11*,nnnnnntFuutFuu1*,11*,2nnnnnnnFFtuutFuub、赫恩格式：（3.73）1*,11*,nnnnnnnBFAFtuutFuu其中，A+B=1。對于（a）有：

16、A=0，B=1； （b）有：A=1/2，B=1/2通用形式：（3.74）1*,11*,nnnnnnnBuAutIuutuIuutBAItBuunn 211tItBuuGnn 211穩定性分析：（3.75）進一步可以將其化為：（3.76）于是：（3.77）1*,11*,nnnnnnnBFAFtuutFuu（1）、當A=0，B=1時，為歐拉后差格式，此時：當滿足 時，格式為穩定的，該格式為條件穩定。（2）、當A=1/2，B=1/2時，為赫恩格式此時：上式表明，無論該式中 取何值， ，故該格式為絕對不穩定格式。tItG 2110111222422tttttG 1t tItG 221142411tG

17、t 1G在上面介紹的時間積分格式均為兩個時間層上的格式，下面我們將簡單介紹包含三個時間層 的時間積分格式。常用的中央差格式，又稱跳背格式，其形式為：（3.78）對于這種格式，可以看出不僅需要一個具有物理意義的初值 ，同樣還需要一個出于計算要求的初值 ，前者稱為物理初值，后者稱為計算初值。3、三個時間層的積分格式11,nnnuuunnntFuu2110u1utIuuuunnnn211tIGG 210122tGIG tItG 211tItG 221nnntuIuu 211據Von-Neumann穩定性判別方法，以上格式可寫成：對應的兩個根為：出現了兩個增幅因子(這表明由于時間積分格式引起了計算解問

18、題)。1t 當 時，格式是穩定的。4 有限差分格式的誤差分析實際工作中，不可能任意縮小步長（由于測站密度的限制，計算量過大等原因），實際計算是在有限的網格下進行的，一定程度的誤差是不可避免的。xt,差分方程可以精確代替微分方程；從而為得到足夠準確的微分方程的近似解提供了保證。0時本節將分析在有限網格區域下差分格式引起的常見誤差。計算解是使用三個時間層積分格式所引起的一個普遍問題。以下以平流方程為例，與其有關的振幅方程為一、計算解問題分析采用三層時間積分格式時，計算解的產生及其特點。uIdtdu tIeutu)0()(一個解若時間微商取中央差格式nnntuIuu 211（3.80）（3.55c）

19、（3.56）可以得到關于增幅因子的二次方程：0122tGIG tItG 211tItG 221其對應的兩個根為：（3.81）1G1lim10Gt1lim20Gt 可見，與 相聯系的波，具有實際的物理意義，稱之為物理解，而與 相聯系的波則不具有任何物理意義，是在計算過程中產生的虛假波型，稱之為計算解。對應的兩個解我們可以表示為：物理解：計算解：1G2G0111uGunn0222uGunn其中a，b為常數。022011ubGuaGunnn而計算所得的數值解為：（3.82）因此，在設計差分格式時，總是使得該差分格式下計算解與物理解振幅的比盡量小，實際工作中通常采用提高網格分辯能力來抑制和減小計算解帶

20、來的影響。以慣性振蕩為例，它反映了旋轉地球上均勻氣壓場中的大氣水平振蕩，其描述方程如下：令 其解為： ，其中 為慣性振蕩頻率。二、時間的截斷誤差（頻率誤差）0, 0fudtdvfvdtduivuWifWdtdWtieWW0f 以下通過數值方法求上述問題的振蕩頻率的差分解（數值解），通過與 比較，說明差分格式所引起的頻率誤差。f 差分方程形為：設其數值解為： （ 為數值解頻率或差分解頻率）。 代入上式，可得：nnntWifWW211tnineeWW0e tfttfttifeeeetitiee1sin1sin2tfttfttWifeWeWeetitiee1000sin1sin2 （a）中央差格式（

21、顯式格式）nnnnWWtifWW112tninieWW0i 121titiiietife tittiftitiiii sin1cos2sin1cos（b）采用梯形格式（隱式格式）差分方程：設其數值解為：（ 為數值解頻率/差分解頻率）代入上式有： 展開之后有：令以上等式兩邊虛部、實部系數相等，則有：2cos1sin1cos2sin2sincos1sin21costfttttfttfttttftiiiiiiii sincos1sincos1sincos1sin22tg 2222tfarctgttfttgii 2cos1sin1cos2sin2sincos1sin21costfttttfttfttt

22、tftiiiiiiii 利用半角公式：故有：即為差分解頻率。可見用時間差商來代替時間微商時會引起頻率誤差。圖3.2時間差分格式的頻率比較tfttfttWifeWeWeetitiee1000sin1sin2 2222tfarctgttfttgii 假設某要素場可以表示為形如 的波動則有：其中k為波數，或稱之為微分解的波數。令 ，則有若用中央差分運算代替微分運算：三、空間的截斷誤差（波數誤差） ikxexuikxikexu xmxxikmxmikxmikmmexxkieexxuusin2121111為差分解波數或計算波數xxksinxikmmeu若用四階精度的差分格式：2211613421mmmm

23、muuuuxxuxxkxk62sinsin8其計算波數為： 為了反映空間差商代替空間微商時引起的波數誤差，定義實際波數計算波數R來表示空間差分格式的精確度。可見，差分格式的精度隨k或 的減小而增大：對于波長較短的波，其產生的波數誤差較大；而波長較長的波，其產生的波數誤差較小。0 xkx1R對于任何波長的波，四階差分格式均比二階差分格式精度高。因此，在實際工作中，縮小格距或采用較高階的差分格式均能提高計算的精確度。空間差分格式四、相速度和群速度誤差0 xActA IkxetutxA,設其形式解為：一維平流方程：Ikcudtdu相速度和群速度的誤差以下以一維平流方程為例，說明差分格式在有限網格下所

24、引起的相速和群速誤差。kc cdkkcddkdcg c為相速，其振動頻率為：其群速度為：可見，對于一維線性平流方程， 為常數，不存在頻散。ccg xIkiietutAuxkxkcIkdtdusinxkxkccsin*xkcdkkcdcgcos*0211xAActAiii對一維平流方程取空間中央差格式，差分方程：假定差分解相速度：差分解群速度：可見：空間差分格式造成了相速度和群速度誤差，從而引起計算頻散。差分格式在波的移動和能量傳播方面均可造成誤差: （1）由于相速度誤差，減慢了平流過程；（2）造成虛假的計算頻散。（1）波長越長，誤差越小；波長越短，誤差就越大；（2）提高網格分辯率，使 取得足夠

25、小，可以提高相 速度，群速度的準確率。0 xkgcc,xxkxkccsin*xkcdkkcdcgcos*差分解相速度：差分解群速度：的誤差也就越小因此：本節討論表明：1、三層時間積分格式存在計算解問題： 計算解對差分解的影響依賴于網格分辨率和波長。2、時間積分格式引起頻率誤差：顯式格式使其頻率明顯增加，振動加快；隱式格式使其頻率明顯減小，振動減慢。3、空間差分格式引起波數誤差： 高階差分格式所引起的波數誤差要比低階格式小； 波長較短的波，誤差尤為嚴重。4、空間差分格式會引起計算頻散： 尤其對于短波，相速度和群速度均會產生很大的誤差。 通常可采用提高網格分辯率的方法減小各種誤差。5 非線性計算穩

26、定性一、計算實例：對一維平流方程：（3.131）（3.132）0 xuutu 022uxtu 0312xuxuutu 10 x第一種形式第二種形式2121114ninininininiuuuuxtuu111116iiiiininiuuuuuxtuu21niniiuuu上述方程可構造不同的差分方程，下面給出其中兩種方案：其中：12（3.133）（3.134）給定兩種不同的初值：1 （3.135） 2（3.136）根據前面有關線性穩定性的討論，其線性穩定性判據為：xiui 2sin0 xiui 2sin5 . 101xtu此時差分格式是滿足線性計算穩定性條件的。假設 004. 0, 1 . 0tx

27、11 . 01 . 0004. 05 . 2maxxtuxtu221iu對兩種差分格式，分別用以上給定的初值進行計算。為了考察二者的穩定性，對每一步的格點總動能：進行了計算，其總動能隨時間的變化如圖所示。 圖3.4 差分格式總動能隨時間的變化曲線（a對應初值1，b對應初值2；實線-方案1，虛線-方案2）以上分析表明： 在滿足線性穩定性的條件下，非線性差分方程仍然會產生不穩定現象，而且不穩定具有突變的特點。 非線性計算不穩定的產生和差分格式有關，初值對計算不穩定具有顯著的影響。以上簡單說明了非線性計算不穩定的產生，下面進一步分析其產生的原因。 有限網格系統能分辨的波的最短波長為 ，對于非線性作用產生的波長小于 的波動，網格系統不能正確地分辨，而把它錯誤地表示成為某種波長大于 的波，從而造成了這種波的誤差混淆誤差。二、混淆誤差-非線性計算不穩定產生的可能原因x2x2x2 我們把長度為 的空間域分為I等分，即 ， 為空間步長，可得I+1個空間格點。函數u在格點上的值可表示為有限級數的形式：（3.137） 由于 ，上式也可寫成： 其中為k波數，而 為網格系統可以分辨的最大波數。混淆誤差的產生：xLxILxx2/02sin2c