第三章_数值计算方案

差分方程及有关的基本概念
若时间和空间都采用中央差分近似，则可以得到如下的差分方 程：
时刻t=（k+1)Δt 在格点x=nΔx 上的值fk+1n ，可以根据t=（k-1)Δt和 t=kΔt时刻的值求出，只要给出初始时刻k=0 的值及边界上的值，问 题可获得解决，(k=1时刻的值可以通过采用时间的前插格式利用初 值计算得到)。上述格式涉及三个时间层，通常称为蛙跃格式。
则称差分格式是稳定的。其中C为与Δt和格距d无关的 参数， εk 是 εk 的某种意义下的范数。
显然一个不稳定的差分格式得到的解是不可能收敛到微 分方程的准确解的，因此设计稳定的差分格式就显得非 常重要。
计算稳定性的分析方法
通常分析差分格式的稳定性有以下几种方 法：
（1）直接证明差分格式的有界性； （2）采用能量法来证明差分格式的稳定性； （3）采用谐波分析的方法。 其中采用谐波分析的方法最为有效和使用方
数值天气预报
第三章 数值计算方案
兰州大学大气科学学院
数值计算方案
通常的大气数值模式的控制方程组是一组复杂的非线性偏微分方 程组，不可能找到一个普遍的解析求解方法，只能采用数值方法 求离散方程的近似解。
常用的数值方法有：
（1）差分方法：采用差商代替微商，使得偏微分方程组变成差分 方程组，可以用代数方法求解；
u(nΔx,kΔt) 和其相容的差分格式的解 unk 之差：
unk − u(nΔx, kΔt)
如果 Δx, Δt → 0 时，在解域满足 max unk − u(nΔx, kΔt) → 0
则称格式是收敛的，即差分格式的解收敛到微分方 程的解。
采用差分方程来近似代替微分方程，这自然会出现误差。 造成数值计算的误差为：
2
if 1 ≥ β ≥ 0 ⇒ G 2 ≤ 1
可以看出 β ≤ 1 是该格式稳定的充分条件。
普遍意义下谐波法分析法
差分方程为：
F k +1 n
=
L( Fnk
)，设初值为：Fn0
=
A0eiμ (nΔx) ,
解的形式为：Fnk
=
Aeiμ (nΔx−μkΔt )
=
A Geiμ (nΔx−μkΔt ) 0
一维线性平流方程的准确解是F(x,t)=F(x-ut,0) ，它表
明 (x,t) 空间任意一点P的值只由初始点Q的值确定，Q
点是通过P点的特征线(t-tp=x-xp)与X轴的交点，对于差 分方程，P 点的解也存在一个依赖区，不难推知它就是
图中AB 线段内的点，直线AP和BP的斜率分别为 r1 = Δt Δx
假定在初始时刻n=n0点有一个误差ε ，其余各点及计算过 程中不再产生误差，考察误差ε将如何传播。不难推得误 差ε满足如下差分方程：
由公式可知，随着时间的增加，误差也增长起来，这种 现象就是计算不稳定现象。因此可以这样来定义计算稳 定性：
差分方程的初始误差 ε0 ，对于任意时间t=kΔt 总有 εk ≤ C ε0
设f
k n
=
Ak eiμxn ,带入差分方程有：Ak+1
= (1− β ) Ak
+ β Ak e−iμΔx
设增幅因子为：G =
Ak +1 Ak
,当
Ak +1 Ak
=
G
≤ 1时格式稳定。
G=(1− β ) + β e−iμΔx ⇒
G 2 =1-2β (1− β )(1− cos μΔx) = 1− 4β (1− β ) sin2 μΔx
为球面函数
Plm (sin ϕ) 为归一化的缔合Legendre多项式
利用(a)的展开方法把所有的预报变量在球面上展开，带入大 气的基本运动方程组，对谱系数进行预报，然后通过逆变换 得到预报场。这种方法的计算量很大。 由于边界比较难于处理，因此很少用来做区域的数 值天气预报和模拟。
用有限项谱展开表示连续函数 线性微分运算可针对基函数直接进行 非线性项采取变换方法（谱模式得以 实现的关键）
和 r2 = − Δt Δx

直线PQ 的斜率等于
r0 =
1 u
。如果
Δt ≤ 1 Δx u
则Q点落入
AB内，反之Q将落到AB 外。如果让 Δt 和 Δx 均趋于0
但始终保持 Δt ≥1 ，Q点的初值与AB中点的初值毫 无关系，也就是Δx说差分方程的解和微分方程的解无关，
不论 Δt 和 Δx 取得如何小，差分解也不会收敛到微分
能保证收敛
稳定性：当 n → ∞ 时数值解如果无穷增大则不稳定
守恒性：差分格式应该尽可能的满足某些物理量的守恒性
质。
色散性：差分格式的解尽管在稳定性的保证下能收敛于源
方程的准确解,但是在准确剧变的位置常常呈现光滑现象,这 种格式称为具有耗散性。
弥散性：在准确剧烈变化的地方还常常呈现微小的”高频振
采用同样的方法可以讨论稳定性问题。
对于三个时间层的格式：Fnk+1 = L(Fnk , Fnk−1)
可以设H
k n
+1
=
Fnk -1，可以化为方程组的形式进行讨论。
离散化有关的概念
相容性：差分系统在某种程度上近似于微分系统
精确性：计算误差问题，截断误差和舍入误差
收敛性：差分方程的解能否收敛到微分方程的解，相容不
设增幅因子为：G =
Ak +1 Ak
,当
Ak +1 Ak
=
G
≤ 1时格式稳定。
对于差分方程组
:
⎧⎪ ⎨ ⎪⎩
F k +1 n
H k +1 n
= =
L1
(
Fnk
,
H
k n
)
L2
(
Fnk
,
H
k n
)
,
解的形式为：⎪⎨⎪⎩⎧FHnknk
= =
A eiμ (nΔx−μkΔt ) 1
A eiμ (nΔx−μkΔt ) 2
e = u − uˆ = u − u + u − uˆ
其中，u 为真实解， uˆ 为数值计算 解， u 为差分方程
的准确解。可以看出造成模式计算误差的两个主要来源 是： （1）截断误差
u − u ，它是由差分方程模拟微分方程的近似中造成
的，这种误差依赖于时间步长和空间格距； （2） 舍入误差
u − uˆ ，它是由计算机的精度造成的。舍入误差在长
用分析稳定性的方法,把改写成
un+1 = 1+ iαωΔt un 1− iβωΔt
得到增幅因子G的表示式
G
=
1+ iαωΔt 1− iβωΔt
=
1−αβ (ωΔt)2 + iωΔt 1+ β 2 (ωΔt)2
z精度的表示：舍去项的阶数
离散化的方法：虽然正方形覆盖计算精度比三角形差些，但其 在程序化方面最容易实现，因此一般的网格覆盖大都采用正方 形网格。设计算区域为二维矩形区域，将它划分为网格：
这样可以用网格点标号(n,m,k) 来表示时空域 (x,y,t) 中点的位 置。由于气象场具有很好的光滑性，因此可以函数f在x点展开为 泰勒级数的形式。以一元函数f(x) 为例有：
荡”,这种格式称为具有色散性
时间积分格式
二时间层积分格式(非迭代格式)
1.格式及其截断误差(R) (1)欧拉(Euler)格式
un+1 = un + ΔtF n , R = O(Δt)
(2)后差格式(隐式格式)
un+1 = un + ΔtF n+1, R = O(Δt)
(3)梯形格式(隐式格式)
u n+1 = u n + Δt ( F n + F n+1 ),
解。这就说明相容并不能保证收敛。
∂u + c ∂u =Biblioteka Baidu0 ∂t ∂x
关于收敛性有一个著名的Lax 等价定理： 对于一个适定的初值问题，如果差分格式是相容的，那么
计算的稳定性是收敛性的充分必要条件。 一般来说，差分格式的收敛性指的是：在解区的所有
网格上，固定nΔt，考察微分方程的初值问题的真解
精度表示：截取的展开项数
截断方式：菱形截断——构造简单 三角截断——各向同性
（经纬向相同分辨率）
（2）有限元法－－－－在气象上很少用
（3）有限差分法－－大气格点模式
对于充分光滑的函数，其导数是当自变量的 增量趋于0 时，因变量的增量与自变量的增 量之商的极限。对于大气系统而言，实际使 用的网格距Δx 与方程所描述运动的特征尺 度相比，是一个很小的量，因此可以采用变 量的增量与自变量的增量之商来近似表示导 数，也即是采用与求导数相反的过程得到相 应的差分方程。
离散化的网格（网格覆盖）： 正方形，正三角形，正六边形，多边形等
三角形精度较好但算法比较复杂，正六边形的精度与正方形的 大体相当，算法比正方形的复杂。因此，通常使用正方形的网 格。 要点：等距（或不等距）节点上的函数值的集合表示连续函数
有限差分代替微分（为了改进差分算子的性能，常对 微分方程先做一定变形）
中尺度模式通常是有限区域模式
对描述大气运动方程的离散化采用以下几种方法
（1）谱方法－－－－大气谱模式
这种方法通常来形成全球或者半球区域的数值模式
J m +J
∑ ∑ φ(λ,ϕ,t) = a2
φlm (t)Ylm (λ,ϕ)
m=−J l= m
（a)
Ylm (λ,ϕ ) = Plm (sin ϕ )eimλ
其中O（Δx）表示余项。引入格点标号，则可以写成：
前差格式 后差格式 中央差分格式
四阶精度差分格式
这种四阶精度的差分格式利用到了n 周围四个网格点的函数 值，增加了计算量，同时也增加了边界条件的处理，因此经 常采用的是中央差分格式。
用差商代替微商所产生的误差称为截断误差。
可以采用如下方法来估计截断误差： 设f 为谐波函数f(x)=Asin(μx), 那么中央差商与微商的比值为：
（2）谱方法：利用适当的基函数（如球谐函数），把解展开成有 限项的线性组合，将对一个变量预测的问题转化为预报展开系数 的问题；
（3）有限元方法：把偏微分方程问题变成相应的泛函极小问题， 以变分原理为基础，又吸收差分方法的思想而发展起来的新方 法。
在这三种方法中，差分方法最为简便，应用最广泛；有限元方法 在大气数值模式的设计中使用较少。本章主要介绍差分方法的基 础知识。有关谱方法将在后面的章节中加以介绍。
利用差分格式得到的解是否可以看作微分方程的近似解呢？或者说
时，差分方程的解是否会逼近于微分方程的解？这涉及
差分解的收敛性。 显然收敛的一个必要条件是
时，差分格式 微分方程，
这称为差分方程与微分方程相一致或者称为相容性。前面给出的几
种差分格式作为微商的近似当然能够保证相容性，但是满足相容性
是否能够保证解的收敛性？
差分方程及有关的基本概念
差分方程：将微商近似用差商代替得到的代数方程称为差分方程。
一维平流方程
针对这样一个微分方程中的每一个微分算子，可以给出若干个差分 格式，从而形成不同的差分方程。对于一个已知的函数而言，在某 个点微分值是唯一确定的，采用了不同的差分格式后，由于不同的 差分格式具有不同的内在性质和与原来微分算子有不同的近似程 度，呈现不同的数值效应。 因此，要使得采用的差分格式得到的数值解能够很好地反映真实的 物理现象，仔细分析各种格式的有效性、可靠性以及长期数值积分 的收敛性显得尤为重要。下面我们给出涉及差分格式的有关基本概 念。
期的数值积分过程中是否增长，这涉及到差分格式的稳定 性问题。
这里利用一维的线性平流方程来讨论。其微分方程 的近似解需要利用差分方程的时间积分才能得到某 一时刻的解。
如果初始时刻或者计算过程中有一个误差，那么这 个误差会不会增长起来？为了简单起见，我们对一 维的线性平流方程的时间微分项采用前插格式，空 间微分项也采用前差格式，则有如下差分方程：
由表可见，对于波长为4Δx以下的波，误差很大，要求相对误差 低于10%，波长须在8Δx以上。
对于二阶导数，可以类似的得到：
对于二元函数的拉普拉斯算子，则有如下的差分格式：
对时间域也可以做类似的离散化。实际的大气物理量 场是在时空域里的，一般我们把时域标号放在上标位 置，而把空域标号放在下标位置以区分时域和空域。
便。
谐波法测试差分格式的有界性
一维线性平流方程及其差分格式为：（u=常数）
∂f + u ∂f = 0方程的准确解是：f (x,t) = F (x − ct) = Aeiμx ∂t ∂x
f k +1 n
−
Δt
f
k n
+u
f
k n
−
fk n−1
0⇒
Δx
f k +1 n
=
(1 −
β
)
fnk
+β
fk n−1
R = O(Δt)2
2
其通用形式写如
u n+1 = u n + Δt(α F n + β F n+1 )
式中α+β=1,当α=1,β=0时为欧拉格式;当α=0,β=1时为后差格式;而当 α=β=1/2时,则是梯形格式.
2,格式分析 根据振动类方程通式可以写成
un+1 = un + iωΔt(αun + β un+1)