第三章 数值计算方案

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