《数值天气预报》中非线性计算不稳定的举例

开篇：润墨网以专业的文秘视角，为您筛选了一篇《数值天气预报》中非线性计算不稳定的举例范文，如需获取更多写作素材，在线客服老师一对一协助。欢迎您的阅读与分享！

摘要：数值天气预报是大气科学的一个重要分支，是一门实用性很强的应用基础学科。数值天气预报本质是用数值方法求解非线性的大气运动方程组。但当采用格点差分来表示微分方程中的非线性项时，易产生非线性计算不稳定现象。本文以一个简单的一维非线性平流方程的数值求解过程为例，给出隐式格式差分方程的基本计算方法并重新演示了非线性计算不稳定现象。由于举例更为简单，可加深学生对显式格式和隐式格式优缺点的理解。

关键词：大气科学；非线性计算不稳定；隐式格式

中图分类号：G642.0 文献标志码：A 文章编号：1674-9324（2013）42-0160-02

一、引言

根据物理定律，比如牛顿第二定律、质量守恒定律、能量守恒定律、气体试验定律等，可以得到支配大气运动的基本方程组。但由于大气运动基本方程组是一组高度非线性的偏微分方程组，很难求得其解析解，人们可通过数值的方法求其近似解（数值解），这就是数值天气预报。数值天气预报是一门实用性很强的应用基础学科[1]。通过差分方法求解大气运动基本方程组时，人们发现数值积分过程中会产生计算不稳定问题，这就需要采用恰当的差分格式或积分格式。在《数值天气预报》课程关于时间积分格式的讲解中[2]，重点介绍了“显示”格式，即差分方程的右端项全部为当前（或和过去）时刻的变量值，通过积数值分可求出方程左端的未来时刻变量值。但对于“隐式”及“半隐式”格式只是简单提及其概念和特点，比如隐式格式为用未来时刻变量值求出未来时刻变量值，它具有计算稳定、但计算复杂的特点。

在《数值天气预报》课程中，需讲解大量的公式推导和讲解，若不能配合简单而又形象的举例和图形，学生尤其是本科生作为授课对象，将很难理解和接受本课程中的相关内容，讲课的效果也将大打折扣。在非线性不稳定计算的举例中，教材中[2]虽给出了采用不同的初值和不同的差分方案对计算稳定性的影响，但并没有清楚地列出其求解过程，因此学生很难了解隐式格式差分的具体求解过程，对教材中列出的“显示”和“隐式”格式的各自优缺点更是难以理解。因此，需要对两种格式的计算过程进行相应的讲解，尤其是隐式格式。此外，教材中[2]对同一个微分方程构造的两个不同的差分方程中，除隐式格式和显式格式的差异外，还存在着对■采用了不同的差分格式，即显式格式采用中央差格式，隐式格式采用前差格式。本课程[2]已清楚的讲解到中央差格式虽具有较高的计算精度，但在时间差分计算时存在计算解的问题，若初值取得不当，则计算解会有较大的振幅。因此，从逻辑上讲，教材中给出的不同的差分方案的影响，实际上不仅仅来源于显式格式和隐式格式的差异，还来源于对时间微分采用不同差分格式的差异。这又加大了学生对显式格式和隐式格式特点的理解难度。

针对上述问题，本文将以简单的一维非线性平流方程为例，给出隐式格式差分方程的具体求解过程，重新探讨非线性计算不稳定现象，目的是使学生更好地了解显式格式和隐式格式差分方程的求解过程，深刻理解两种格式各自的优缺点。

二、非线性计算不稳定的计算实例

以大气科学中极具代表性的一维平流方程为例：■+u■=0，0≤x

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

或■+■（u■+■），0≤x

以上两式与教材[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

三、隐式格式的求解

显式差分方程（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.